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ABSTRACT 


'  Until  very  recently  the  Transit  satellite  positioning  system 

was  the  only  accurate  navigation  system  available  on  a  (practically) 
global  basis.  Therefore  Transit  is  a  commonly  used  navigation  aid  and 
4  was  included  in  the  DREO  developed  Marine  Integrated  Navigation 

System  (MINS-B  II)  and  Primary  Land  Arctic  Navigation  System  (PLANS). 
MINS  and  PLANS  are  both  Kalman  filter  based  multi-sensor  optimally 
integrated  navigation  systems  and  as  such  require  detailed  stochastic 
and  deterministic  error  models  for  all  of  their  sensors.  This  paper 
provides  a  detailed  geometric  derivation  of  the  deterministic  error 
model  relating  the  position  error  of  a  Transit  satellite  Doppler 
position  fix  (or  similarly  a  SARSAT  fix)  to  the  error  in  the  velocity 
and  height  that  is  fed  into  the  Transit  receiver  (or  assumed  for  the 
SARSAT  transponder).  This  error  model  is  presented  in  the  form  of 
explicit  functions  of  variables  which  are  normally  provided  by  a 
Transit  receiver  with  each  position  fix,  namely  the  satellite  maximum 
elevation  angle,  the  satellite  direction  of  travel  and  the  receiver 
latitude.  Sample  plots  of  these  sensitivity  functions  are  presented, 
along  with  experimental  results  for  verification.  This  model  is  used 
in  MINS  and  PLANS  to  form  the  Kalman  filter  measurement  error 
sensitivity  matrices. 


R&SUH& 

Jusqu'a  tres  r^cemment,  le  systeme  de  posit ionnement  par 
satellite  Transit  6tait  le  seul  systeme  de  navigation  precis 
disponible  pratiquement  partout.  Par  consequent,  Transit  est  un  aide 
a  la  navigation  utilise  communement  et  a  et£  inclus  dans  le  Systeme 
de  Navigation  Integre  pour  la  Marine  (SNIM-B  II)  et  pour  le  Systeme 
de  Navigation  Primaire  pour  la  Terre  Arctique  (SNPTA)  developpes  par 
le  CRDO.  SNIM  et  SNPTA  sont  tous  deux  des  systemes  de  navigation 
integres  da  fa?on  optimum  utilisant  de  multiples  d£tecteurs  bases  sur 
le  filtre  Kalman,  et  commes  tels,  requierent  des  modeles 
stochastiques  et  deterministiques  detailles  par  les  erreurs  pour  tous 
les  detecteurs.  Cet  article  fournie  une  derivation  geometrique 
detaillee  d'un  modele  pour  l'erreur  deterministique  reliant  l'erreur 
sur  la  position  d'un  satellite  Transit  par  un  fix  Doppler  (ou  de 
faqon  similaire  un  fix  SARSAT)  a  l'erreur  sur  la  vitesse  et  la 
hauteur  qui  est  communique  au  receveur  Transit  (ou  assumee  pour  le 
transpondeur  SARSAT).  Ce  modele  pour  l'erreur  est  presente  sous  la 
forme  de  fonctions  explicites  de  variables  qui  sont  normalement 
fournis  par  le  receveur  Transit  avec  chaque  fix  sur  la  position,  e.g. 
l'angle  d'614vation  maximum  du  satellite,  la  direction  de  la 
trajectoire  du  satellite  et  la  latitude  de  receveur.  Des  graphiques 
Achantillons  de  ces  fonctions  de  sensibility  sont  presentes,  avec  des 
r£sultats  exp^rimentaux  pour  v6rif ication.  Ce  modele  est  utilise  pour 
SNIM  et  SNPTA  pour  former  des  mesures  du  filtre  Kalman  pour  des 
matrices  de  sensibility  des  erreurs. 
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EXECUTIVE  SUMMARY 


Transit  is  a  satellite  based  navigation  system,  which  was 
originally  developed  for  the  U.S.  Navy  Polaris  submarine  fleet  by 

*  the  Applied  Physics  Lab  of  John  Hopkins  University.  The  system 
has  been  operational  since  1964,  and  was  released  for  public  use 
in  1967. 

The  Transit  Satellite  positioning  system  basically 
consists  of  5  or  6  satellites  in  circular  polar  orbits,  with  an 
orbital  height  of  about  1,075  km,  and  a  107  minute  period.  These 
satellites  transmit  continuously  at  two  very  stable  frequencies 
(150  MHz  and  400  MHz).  An  earthbound  receiver  can  obtain  a 

position  fix  whenever  a  satellite  passes  overhead  by  measuring 
the  Doppler  frequency  shifts  due  to  the  relative  motion  as  the 
satellite  rises  towards  the  receiver,  passes  over  and  sets  away 
from  the  receiver.  The  transmitted  signal  is  modulated  with  a 
data  message  which  contains  ephemerides  and  timing  marks.  From 
this  message  the  receiver  can  accurately  calculate  the  absolute 
position  and  velocity  of  the  satellite  throughout  the  pass.  From 
this  known  satellite  position  and  velocity  profile,  and  the 
Doppler  derived  relative  velocity  profile,  the  receiver  can 

»  calculate  its  own  position. 

In  order  to  calculate  its  position  however,  the  receiver 

*  must  be  able  to  remove  the  effect  of  its  own  velocity  from  the 

Doppler  measurement.  Thus  the  receiver  must  either  be  stationary 

during  the  satellite  pass,  as  is  the  case  with  survey 

instruments,  or  the  receiver  motion  (velocity)  relative  to  the 
earth  must  be  known  during  the  pass.  Generally  then,  the  receiver 
must  be  continuously  fed  the  receiver  velocity.  Any  error  in  this 
velocity  input  will  lead  to  an  error  in  the  position  fix  that  the 
receiver  produces.  The  exact  relationship  between  this  velocity 
input  error  and  position  output  error  is  the  primary  concern  of 
this  paper.  This  relationship  is  needed  primarily  for  use  in 
Kalman  filter  based  integrated  navigation  systems,  where  velocity 
and  position  errors  are  estimated.  DRE0  has  developed  two  such 
multi-sensor  systems:  MINS  (Marine  Integrated  Navigation  System) 
and  PLANS  (Primary  Land  Arctic  Navigation  System),  which  use 
Transit  as  one  of  their  sensors. 

There  is  another  important  source  of  error  that  can  be 
modelled  in  a  similar  manner,  namely  the  position  error  that  is 
caused  by  the  error  in  the  height  above  sea  level  that  is 
,  provided  to  the  receiver.  Although  it  is  possible  to  use  the 

Doppler  measurements  to  solve  for  the  three  dimensional  position 
of  the  receiver  (four  dimensional  if  time  is  included),  the 
accuracy  is  very  poor  in  the  direction  that  is  orthogonal  to  the 

*  satellite  velocity  vector  and  to  the  receiver-satellite  range 
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vector.  The  "error  ellipsoid"  is  therefore  very  prolate,  being 
elongated  in  this  one  direction  and  narrow  in  the  others.  This 
would  generally  cause  an  uncertainty  in  all  three  spatial 
dimensions.  Fortunately  it  is  possible  to  solve  this  problem  by 
providing  the  receiver  with  height  above  a  reference  ellipsoid 
("sea  level").  This  height  surface  will  intersect  the  error 
ellipsoid  in  a  small,  nearly  circular  error  ellipse,  (except  in 
the  case  when  the  major  axis  of  the  ellipsoid  is  near  the 
horizontal  plane  in  which  case  the  satellite  maximum  elevation 
angle  is  close  to  90  degrees  and  the  receiver  will  be  unable  to 
solve  for  position).  Therefore  most  Transit  receivers  require  the 
user  to  supply  a  height  input  (a  constant  for  marine 
applications)  and  consequently  any  error  in  this  height  will 
produce  a  position  error. 

This  report  provides  a  complete  derivation  of  the 
explicit  linear  relationship  between  the  velocity  and  height 
error  input  and  the  position  error  output  of  a  Transit  receiver. 
The  derivation  is  geometric  in  nature,  and  provides  some  useful 
insights  along  the  way.  A  strictly  analytical  approach  was  not 
used  because  it  seemed  too  intractable,  and  without  intermediate 
verifiable  steps.  Therefore  the  geometric  approach  was  taken, 
which  although  more  intuitively  complicated,  had  the  great 
benefit  of  having  intermediate  results  which  had  clear  physical 
significance  and  could  be  verified  to  some  extent.  Thus  through  a 
combination  of  elementary  spherical  geometry,  planar  geometry  and 
differential  techniques,  the  required  exact  linear  relations  were 
found  in  the  form  of  a  2x3  matrix,  whose  components  are  explicit 
functions  of  4  parameters  which  are  normally  provided  by  the 
Transit  receiver  with  its  position  fix  (the  maximum  satellite 
elevation  angle,  the  receiver  latitude,  the  satellite  north-south 
directions  of  travel  and  east-west  direction  from  the  receiver). 

An  algorithm  is  given  to  efficiently  calculate  this  error 
relationship.  Numerical  results  are  presented  to  illustrate  the 
functional  dependence  of  the  sensitivity  matrix  elements.  One 
significant  new  result  shown  is  the  substantial  change  in  these 
functions  at  high  latitudes.  The  qualitative  error  models  shown 
in  the  literature  are  therefore  seen  to  be  invalid,  even 
qualitatively,  at  high  latitudes.  Substantial  experimental 
results  are  presented  which  verify  the  high  fidelity  of  the  model 
for  all  satellite  directions  of  travel,  at  all  satellite 
elevation  angles,  at  mid-latitudes. 

One  of  the  intermediate  results  found  in  this  way  was  a 
lower  bound  which  can  be  placed  on  the  maximum  elevation  angle  of 
a  satellite  pass,  as  a  function  of  receiver  latitude.  This 
provides  a  consistency  check  on  the  elevation  angle  provided  as 
part  of  the  Transit  fix  data  set.  Upper  and  lower  bounds  were 
also  found  for  the  range  to  the  satellite  subpoint,  which  can  be 
used  to  justify  linearization. 
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1.0  INTRODUCTION 


Transit  is  a  satellite  based  navigation  system,  which  was 
originally  developed  for  the  U.S.  Navy  Polaris  submarine  fleet  by 
the  Applied  Physics  Lab  of  John  Hopkins  University.  The  system 
has  been  operational  since  1964,  and  was  released  for  public  use 
in  1967.  Reference  [11  gives  a  more  complete  description. 

The  Transit  Satellite  positioning  system  basically 
consists  of  5  or  6  satellites  in  circular  polar  orbits,  with  an 
orbital  height  of  about  1,075  km,  and  a  107  minute  period.  These 
satellites  transmit  continuously  at  two  very  stable  frequencies 
(150  MHz  and  400  MHz).  An  earthbound  receiver  can  obtain  a 
position  fix  whenever  a  satellite  passes  overhead  by  measuring 
the  Doppler  frequency  shifts  due  to  the  relative  motion  as  the 
satellite  rises  towards  the  receiver,  passes  over  and  sets  away 
from  the  receiver.  The  transmitted  signal  is  modulated  with  a 
data  message  which  contains  ephemerides  and  timing  marks.  From 
this  message  the  receiver  can  accurately  calculate  the  absolute 
position  and  velocity  of  the  satellite  throughout  the  pass.  From 
this  known  satellite  position  and  velocity  profile,  and  the 
Doppler  derived  relative  velocity  profile,  the  receiver  can 
calculate  its  own  position.  Section  4.0  below  gives  more  details. 

In  order  to  calculate  its  position  however,  the  receiver 
must  be  able  to  remove  the  effect  of  its  own  velocity  from  the 
Doppler  measurement.  Thus  the  receiver  must  either  be  stationary 
during  the  satellite  pass,  as  is  the  case  with  survey 
instruments,  or  the  receiver  motion  (velocity)  relative  to  the 
earth  must  be  known  during  the  pass.  Generally  then,  the  receiver 
must  be  continuously  fed  the  receiver  velocity.  Any  error  in  this 
velocity  input  will  lead  to  an  error  in  the  position  fix  that  the 
receiver  produces.  The  exact  relationship  between  this  velocity 
input  error  and  position  output  error  is  the  primary  concern  of 
this  paper,  and  is  derived  in  Sections  6.1  to  6.4  below. 

There  is  another  important  source  of  error  that  can  be 
modelled  in  a  similar  manner,  namely  the  position  error  that  is 
caused  by  the  error  in  the  height  above  sea  level  that  is 
provided  to  the  receiver.  Although  it  is  possible  to  use  the 
Doppler  measurements  to  solve  for  the  three  dimensional  position 
of  the  receiver  (four  dimensional  if  time  is  included),  the 
accuracy  is  very  poor  ii:  the  direction  that  is  orthogonal  to  the 
satellite  velocity  vector  and  to  the  receiver-satellite  range 
vector.  The  "error  ellipsoid"  is  therefore  very  prolate,  being 
elongated  in  this  one  direction  and  narrow  in  the  others.  This 
would  generally  cause  an  uncertainty  in  all  three  spatial 
dimensions.  Fortunately  it  is  possible  to  solve  this  problem  by 
providing  the  receiver  with  height  alove  a  reference  ellipsoid 
("sea  level").  This  height  surface  will  intersect  the  error 
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ellipsoid  in  a  small,  nearly  circular  error  ellipse,  (except  in 
the  case  when  the  major  axis  of  the  ellipsoid  is  near  the 
horizontal  plane  in  which  case  the  satellite  maximum  elevation 
angle  is  close  to  90  degrees  and  the  receiver  will  be  unable  to 
solve  for  position).  Therefore  most  Transit  receivers  require  the 
user  to  supply  a  height  input  (a  constant  for  marine 
applications)  and  consequently  any  error  in  this  height  will 
produce  a  position  error. 

This  report  describes  the  linearized  relationship  between 
the  velocity  and  height  error  input  and  the  position  error 
output.  This  relationship  is  expressed  in  the  form  of  a 
sensitivity  matrix  H,  where: 


■s 

north  position  error 

/  > 
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where  •  difference  "A"  is  defined  to  be  the  true  value  minus  the 
estimate.  For  example: 


AN  ■  N  -  N 


As  shall  be  shown,  the  components  of  this  sensitivity 
matrix  H  are  complicated  functions  of  the  satellite  maximum 
elevation  angle,  the  satellite  direction  of  travel  (north  to 
south  or  vice  versa),  the  receiver  latitude,  and  whether  the 
satellite  subpoint  is  east  or  west  of  the  receiver  at  maximum 
elevation.  These  parameters  are  defined  in  the  next  section. 

The  relationship  defined  by  equation  (2)  is  not  only 
useful  in  predicting  accuracies,  but  forms  an  important  part  of 
the  measurement  matrix  for  a  Kalman  filter.  The  measurement 
equation  of  a  Kalman  filter  has  the  form: 
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H(X)  +  V 


(3) 


which  relates  the  state  vector  X,  of  quantities  that  the  filter 
is  trying  to  estimate,  to  the  measurement  vector  Z  of  quantities 
fed  into  the  filter,  and  the  measurement  noise  vector  V.  In  a 
linear  filter  the  vector-valued  function  H(X)  is  simply  the 
product  of  a  "measurement  matrix"  H  with  the  state  vector  X,  and 
the  measurement  equation  becomes: 


Z  =  H  X  +  V  (4) 


At  DREO  we  have  developed  two  Kalman  filter  based 
optimally  integrated  navigation  systems  which  utilize  Transit 
measurements:  a  Marine  Integrated  Navigation  System  (MINS), 
described  in  reference  [2],  and  a  Primary  Land  Arctic  Navigation 
System  (PLANS),  described  in  reference  [3].  It  was  for  these  two 
systems  that  this  detailed  error  model  was  developed.  In  1983  the 
author  first  developed  an  implicit  error  model,  for  use  in  MINS. 
The  numerical  results  were  presented  at  reference  [4]  and  are 
essentially  the  same  as  the  results  shown  in  Section  9.0  below. 
The  explicit  error  model  derived  in  this  report  was  developed  in 
1989,  for  the  PLANS  ADM.  This  explicit  model  leads  to  a  much  more 
efficient  algorithm,  and  greater  confidence  in  the  results. 

This  error  model  is  particularly  important  for  the  PLANS 
application,  where  the  relationship  given  by  equation  (2)  is  in 
effect  partially  inverted  to  solve  for  an  otherwise  unbounded 
velocity  error.  This  velocity  error  arises  in  high  latitude 
operations  between  the  geographic  and  geomagnetic  poles  where  an 
accurate  heading  cannot  be  reliably  obtained,  either 
magnetically  or  with  a  gyrocompass. 

A  geometric  approach  is  taken  in  deriving  the  H  matrix  of 
equation  (3),  and  some  preliminary  geometric  quantities  must  be 
defined  first. 
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2.0  DEFINITIONS 


Table  1  lists  the  relevant  variables  and  parameters  while 
Figure  1  illustrates  the  primary  quantities  of  interest.  Note 
that  vectors  are  underlined.  The  magnitude  of  a  vector  is 
represented  by  the  vector  symbol  without  the  underline. 


Table  1:  Summary  of  Relevant  Variables. 


PRIMARY  VARIABLES:  (whose  relationships  are  sought) 

Av 

receiver  velocity  error  (true  -  estimate) 

Aa 

heading  of  receiver  velocity  error  (not  heading  error) 

Ah 

Receiver  height  error 

Ar 

receiver  position  error  (AN,AE) 

AN 

North  position  error  of  receiver  (due  to  Av,  Aa  &  Ah) 

AE 

East  position  error  of  receiver  (due  to  Av,  Aa  &  Ah) 

SECONDARY:  (available  directly  from  the  receiver  if  needed) 

r 

Receiver  position  vector  (latitude  X,  longitude  L) 

V 

Receiver  speed 

a 

Receiver  heading 

h 

Receiver  height 

a 

Satellite  maximum  elevation  angle 

NS 

satellite  north-to-south  direction  of  travel  (at  PCA) 

(  »  1  if  towards  the  south,  -1  if  towards  the  north  ) 

satellite  east-to-west  direction  from  receiver  (at  PCA) 

(  »  1  if  towards  the  east,  -1  if  towards  the  west  ) 

INTERMEDIATE:  (derived  from  secondary  variables  if  needed) 

s 

Satellite  position  vector 

2 

Satellite  subpoint  position  vector  (Xp,  Lp) 

p 

slant  range  from  receiver  to  satellite 

e 

subpoint/  receiver  separation  angle  (at  earth  centre) 

* 

bearing  (from  north)  to  subpoint  from  receiver 

♦ 

bearing  (from  north)  to  receiver  from  subpoint 

0 

heading  (from  north)  of  subpoint  track  at  subpoint 

Y 

angle  between  slant  range  and  satellite  velocity 

n 

angle  between  slant  range  and  receiver  velocity  error 

V 

angle  from  receiver  velocity  error  to  subpoint  at  r 

V 

satellite  velocity  vector 

Note:  all  these  quantities  are  defined  in  an  earth  fixed  frame. 
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Figure  1.  Primary  Variables  of  a  Transit  Satellite  Pass 
(in  earth  fixed  coordinates). 


Earth  Centre 


SI  -satellite  position  at  time  tl 
Y1  -cone  angle  at  time  tl 
PI  -slant  range  at  time  tl 
r  -receiver  position 
0  -satellite  elevation  angle 
V  -bearing  angle  to  subpoint 
§  -bearing  angle  from  subpoint 
R  -earth  radius 
H  -satellite  orbital  height 
0  -subpoint-receiver  separation  angle 


Table  2:  Relevant  Constants. 


R 

earth  radius 

=  6,370,000 

metres 

H 

satellite  orbital  height 

=  1,075,000 

metres 

Vs 

satellite  "speed" 

=  7,290 

m/s 

<•> 

earth  rate 

=  0.0000729211585  radians/sec. 

This  satellite  "speed"  Vs  is  constant  in  an  earth  centred 
non-rotating  frame  (not  earth  fixed).  When  earth  rate  at  orbital 
height  is  added  to  this,  the  correct  earth  fixed  speed  is  (since 
these  velocity  components  are  orthogonal): 


V  =  /  Vs2  +  (R+H) 2to2cos2 Xp 


(5) 


-  /  7, 290 2  +  543 2 cos 1 Xp  (6) 

thus 

7,290  m/s  <  V  <  7,310  m/s  (7) 


All  quantities  (except  Vs  as  just  explained)  are  defined 
in  an  earth  fixed  reference  frame.  Figure  1  shows  a  receiver  at 
point  r,  moving  with  speed  v  and  heading  a,  a  satellite  moving 
south  to  north  in  a  polar  orbit.  Two  satellite  positions  are 
shown  si  and  s2,  for  two  points  in  time  tl  and  t2.  The  distance 
between  the  receiver  and  satellite  is  pi  and  p2,  and  is  called 
the  slant  range.  Since  the  receiver  velocity  is  assumed  to  be 
much  smaller  than  the  satellite  velocity  (as  it  would  be  in  the 
marine  and  land  vehicle  applications  of  interest  here),  the 
change  in  receiver  position  between  tl  and  t2  is  not  visible 
here  (a  Transit  satellite  pass  lasts  less  than  20  minutes,  so 
that  a  land  vehicle  could  not  move  more  than  about  30  to  40 
kilometers,  and  in  most  cases  would  move  much  less). 

The  satellite  subpoint  is  the  point  £  on  the  surface  of 
the  earth  directly  below  the  satellite.  Every  satellite  pass  has 
a  "point  of  closest  approach"  (PCA),  which  is  the  point  (in  time 
and  space)  at  which  the  slant  range  is  at  its  minimum.  The  path 
of  the  subpoint  on  the  earth's  surface  is  called  the  subpath. 
Since  the  satellite's  orbit  is  fixed  in  an  earth-centred 
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"inertial  space"  (ie.  non-rotating)  the  subpath  does  not  follow  a 
meridian  of  longitude.  The  meridians  are  shown  for  the  receiver 
and  the  subpoint,  in  order  to  define  the  bearing  from  the 
receiver  to  the  subpoint,  the  bearing  $  from  the  subpoint  to  the 
receiver,  and  the  heading  of  the  subpath  |3  (these  are  all 
clockwise  angles  from  north).  Also  illustrated  is  the  separation 
angle  0,  at  the  centre  of  the  earth  between  the  receiver  and  the 
subpoint. 

It  will  be  seen  that  a  very  important  variable  is  the 
"elevation  angle”  a,  which  is  the  angle  at  the  receiver  from  the 
slant  range  vector  to  the  locally  horizontal  plane.  As  shall  be 
seen,  the  maximum  elevation  angle  corresponds  to  the  PCA.  Another 
important  variable  is  the  "satellite  direction  of  travel"  which 
is  really  two  logical  variables  simply  indicating  whether  the 
satellite  is  moving  north  to  south  or  vice  versa  (NS),  and 
whether  the  PCA  is  east  or  west  of  the  receiver's  meridian  (EV). 
If  the  PCA  is  on  (or  very  close  to)  the  receiver's  meridian  then 
the  geometry  does  not  allow  a  position  fix  to  be  made. 

Some  symbolic  conventions  will  be  used:  in  general  the 
error  in  a  quantity  x  is  represented  by  using  the  prefix  A,  as  in 
Ax,  a  differential  error  is  represented  by  the  prefix  A,  as  in 
&x,  and  an  infinitesimal  tine  difference  is  represented  by  the 
prefix  d,  as  in  dx. 
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3.0  ASSUMPTIONS 


In  this  mathematical  derivation,  all  important 
assumptions  shall  be  explicitly  stated,  and  any  approximations 
clearly  pointed  out. 


In  order  to  apply  spherical  trigonometry  it  shall  be 
assumed  that  the  satellites  are  all  in  circular  polar  orbits  and 
that  the  earth  is  spherical.  The  earth  is  actually  slightly 
oblate,  with  the  best  fitting  ellipsoid  (World  Geodetic  Survey 
1984)  having  semi  major  axis  A  and  semi  minor  axis  B  of  lengths 


A  =  6,378,137  metres 
B  =  6,356,752  metres 


(8) 


respectively,  therefore  having  a  small  eccentricity  of 

f  2  2]1/2 

e  =  l  A  -  BJ  =  0.08  (9) 

A 


This  is  close  enough  to  being  spherical  for  the  purpose  of 
determining  the  desired  error  relationship,  which  is  essentially 
a  local  phenomenon.  The  eccentricity  of  a  typical  Transit 
satellite  orbit  is  even  smaller,  at  about  0.002.  The  deviation 
from  a  polar  orbit  is  usually  expressed  as  the  orbital 
inclination  angle,  which  is  the  angle  between  the  orbital  plane 
and  the  equatorial  plane.  For  Transit  this  is  typically  well 
within  half  a  degree  of  polar  (90  degrees). 

To  obtain  the  linear  relationship  of  equation  (1),  using 
differential  techniques,  it  shall  be  assumed  that  the  Transit 
position  fix  is  made  using  relative  velocity  measurements 
(Doppler)  in  a  neighbourhood  of  the  point  of  closest  approach 
(PCA).  The  PCA  is  a  point  in  time  that  typically  occurs  near  the 
centre  of  the  satellite  pass,  when  the  distance  between  the 
satellite  and  the  receiver  is  a  minimum  (hence  the  satellite 
velocity  in  the  receiver  reference  frame  is  orthogonal  to  the 
receiver-satellite  range  vector).  Receivers  actually  use  Doppler 
measurements  over  the  duration  of  the  satellite  pass,  which  is 
from  10  to  16  minutes,  typically  forming  one  Doppler  count 
measurement  every  23  seconds.  Although  this  is  the  most  serious 
simplifying  assumption,  it  is  justified  on  the  grounds  that  the 
satellite  pass  is  approximately  symmetric  about  the  PCA  and 
therefore  the  effects  of  the  different  geometry  on  the  two  sides 
of  the  PCA  should  approximately  cancel. 
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This  appears  to  carry  the  implicit  assumption  that  the 
velocity  and  height  errors  are  constant  over  the  satellite  pass, 
but  this  is  really  not  a  restriction  if  we  use  the  average  error 
values  in  (1).  The  linearity  of  (1)  permits  averaging  on  both 
sides. 


It  is  also  necessary  to  assume  that  the  receiver  velocity 
error  is  small  compared  to  the  satellite  velocity  (  v  «  V  ) 
which  in  light  of  (7)  is  not  a  serious  restriction  for  the  usual 
terrestrial  and  marine  applications.  This  is  to  ensure  that  the 
error  relationship  of  equation  (3)  is  linear,  as  in  (1). 

There  are  many  equivalent  ways  to  calculate  receiver 
position,  given  the  satellite  position  and  Doppler  measurements 
over  an  interval.  In  this  report  a  largely  geometric  rather  than 
analytic  solution  shall  be  used,  where  the  various  sources  of 
error  and  the  intermediate  quantities  that  must  be  calculated, 
have  clear  physical  interpretations  that  are  easily  understood. 
This  also  facilitates  verification. 
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4.0  TRANSIT  POSITIONING 


4.1  Stationary  Receiver 


For  a  stationary  receiver  r,  in  an  earth  fixed  frame, 
Figure  2  illustrates  geometrically  how  a  Transit  position  fix  can 
be  obtained. 

A  Transit  Doppler  measurement  gives  the  rate  of  change  of 
the  range  from  the  receiver  to  the  satellite.  From  the 
transmitted  ephemerides  data  the  receiver  knows  the  satellite 
position  s,  and  velocity  V,  in  an  earth  fixed  reference  frame.  If 
the  receiver  position  r  is  not  moving  in  this  frame,  then  the 
measured  relative  velocity  (the  rate  of  change  of  the  slant 
range)  dp/dt  is  due  entirely  to  the  satellite  velocity  V,  which 
must  therefore  be  at  an  angle  y  to  the  range  line,  as  shown  in 
Figure  2,  where 


p  *  -V  cosy 


(10) 


(where  a  dot  ’  indicates  d/dt).  Equation  (10)  can  then  be  solved 
for  y.  Then  r  must  lie  on  the  cone  defined  by  the  angle  y  and  the 
satellite  velocity  vector  V.  Two  Doppler  measurements  therefore 
define  two  such  cones,  at  angles  yl  and  y2  about  the  respective 
satellite  velocity  vectors  VI  and  V2.  The  intersection  of  these 
two  cones  is  a  curve  C,  which  becomes  circular  about  the 
satellite  in  the  limit  as  the  time  between  measurements  becomes 
small  (so  that  VI  approaches  colinearity  with  V2). 

Assuming  that  r  is  on  the  earth's  surface,  or  a  surface 
of  known  height,  then  r  must  lie  on  the  intersection  of  this 
surface  and  the  curve  C. 

This  intersection  (of  a  circle  and  a  sphere)  will 
generally  consist  of  two  points,  (provided  that  the  satellite 
does  not  pass  directly  over  the  receiver,  in  which  case  the 
circle  will  be  tangent  to  the  sphere).  The  ambiguity  can  be 
resolved  by  a  third  Doppler  measurement,  because  the  satellite 
path  is  not  coplanar  in  the  earth  fixed  reference  frame.  Since  we 
are  interested  here  only  in  the  effect  of  receiver  velocity  error 
on  the  position  fix,  it  shall  be  assumed  that  this  ambiguity  has 
been  resolved. 

In  fact  these  multiple  Doppler  measurements  can  in 
principle  be  used  to  solve  for  receiver  height  above  geoid  as 
well  by  intersecting  different  curves  C.  However  the  vertical 
geometric  dilution  of  precision  (VDOP)  is  generally  such  that  the 
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height  obtained  from  a  single  pass  would  be  very  inaccurate  since 
the  curves  C  intersect  at  a  sharp  angle.  This  would  degrade  the 
horizontal  solution  since  the  curves  C  are  not  necessarily 
vertical.  Therefore  it  is  normal  for  Transit  receivers  to  be 
given  this  height  information  from  some  external  source. 
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4.2  Moving  Receiver 


If  the  receiver  is  moving  in  this  earth  fixed  frame  then 
there  are  several  complications.  First  the  measured  Doppler  shift 
must  be  compensated  to  remove  the  effect  of  receiver  velocity 
before  calculating  the  angles  yl  and  y2.  This  is  accomplished  by 
subtracting  the  component  of  the  receiver's  velocity  in  the 
direction  of  the  slant  range  vf  from  the  measured  velocity  dp/dt. 
Equation  (10)  then  becomes: 


-V  cosy  =  p  -  vr  (11) 

The  receiver  velocity  component  vf  is  shown  in  Figure  3  to  be: 

vf  *  v  cos(^  -  a)cos<r  (12) 

where  v  and  a  are  the  receiver  speed  and  heading,  and  ^  and  a  are 
as  defined  in  Section  2.0  and  shown  in  Figure  1  (the  bearing  to 
the  subpoint  and  the  elevation  angle). 

Secondly  the  receiver's  position  will  have  changed 
between  times  tl  and  t2,  and  will  therefore  not  be  on  the 
intersection  of  the  cones.  This  can  also  be  resolved  by  applying 
the  appropriate  translation  to  the  first  cone,  as  shown  in  Figure 
4.  The  appropriate  translation  is  of  course  the  position  change 
vector  of  the  receiver,  which  for  a  small  time  interval  is  vAt. 
This  therefore  also  requires  knowledge  of  the  receiver  velocity 
v,  and  in  particular  its  component  in  the  direction  orthogonal  to 
the  slant  range. 
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5.0  EFFECT  OF  ANTENNA  HEIGHT  ERROR 

The  position  fix,  as  explained  in  Section  4.1,  is  found 
by  intersecting  the  Doppler  curve  C  (defined  by  the  intersection 
of  two  Doppler  cones  as  seen  in  Figure  2)  with  a  surface  of  known 
height  h  above  the  earth's  surface.  The  effect  of  an  error  in 
this  height  Ah,  simply  has  the  effect  of  moving  the  position  fix 
along  the  curve  C.  If  the  height  error  is  small  compared  to  the 
satellite  range,  then  this  curve  can  be  approximated  by  its 
tangent  vector  at  the  receiver  location,  as  shown  in  Figure  5, 
where  the  correct  position  is  r  and  the  error  in  the  estimated 
position  is  Ar. 

The  effect  of  Ah  is  then  linear,  in  the  direction  of  this 
tangent  vector.  Thus  this  simple  geometric  effect  can  be 
described  in  terms  of  the  "bearing"  and  "elevation  angle"  of  this 
tangent  vector. 

Assuming  a  circular  orbit,  the  satellite  velocity  V  is 
perpendicular  to  the  vertical  line  sp  from  the  satellite  to  its 
subpoint.  At  the  PCA  the  satellite  velocity  V  is  also 
perpendicular  to  the  slant  range  line  sr  from  the  satellite  to 
the  receiver  (by  definition  of  PCA).  Therefore  at  the  PCA,  V  is 
perpendicular  to  the  plane  srp  containing  the  satellite,  its 
subpoint  and  the  receiver. 

Since  the  Doppler  curve  C  is  perpendicular  to  V  at  PCA, 
and  passes  through  r,  it  must  also  lie  in  the  plane  srp. 
Therefore,  as  shown  in  Figure  5  (where  p,  C  and  U  are  coplanar) 
the  angle  £  between  the  Doppler  curve  C  and  the  vertical,  is 
equal  to  the  elevation  angle  a  (since  both  are  complementary  to 
the  angle  between  the  slant  range  and  the  vertical).  Also  the 
bearing  angle  from  Ar  to  north  is  (*<>-n)  (for  a  positive  Ah,  or 
just  ik  for  a  negative  Ah),  where  is  the  bearing  to  the  subpoint 
path. 

Therefore  the  magnitude  of  the  position  error  Ar  caused 
by  a  height  error  Ah,  satisfies: 


tana  =  tan£  =  Ar/Ah  (13) 


(where  this  Ar  is  defined  to  be  positive  away  from  the  subpoint, 
and  negative  towards  the  subpoint,  giving  it  the  same  sign  as 
Ah).  Thus 


Ar  =  Ah  tana 


(14) 
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Taking  components  o£  Ar  in  the  north  and  east  direction,  we  have: 

AN  ■  Ah  cos(t-n)  tana 

=  -Ah  cosij/  tana  (15) 

AE  a  Ah  sin(ij/-n)  tana 

*  -Ah  sinij/  tana  (16) 

This  gives  us  the  first  column  of  the  H  matrix  of  equation  (2): 


'AN' 

.AE, 


-costj/  tana' 
.-sinij/  tana, 


(17) 


(in  the  absence  of  velocity  errors).  The  satellite  elevation 
angle  a  is  provided  by  the  receiver.  The  bearing  angle  to  the 
subpoint  can  be  expressed  in  terms  of  the  known  quantities 
(X,  a,NS,EW)  as  shown  in  the  appendices.  Appendix  A  shows  how  *(/ 
can  be  expressed  in  terms  of  the  intermediate  variables  (♦, Xp) 
and  appendices  A,  B  and  C  show  how  these  intermediate  variables 
can  be  expressed  in  terms  of  the  known  quantities. 
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Figure  5.  Effect  of  Height  Error. 
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6.0  EFFECT  OF  VELOCITY  ERROR 


The  error  in  the  Transit  position  fix,  due  to  the  error 
in  the  receiver's  estimate  of  its  velocity,  is  decomposed  into 
four  separate  and  independent  errors.  As  was  explained  in  Section 
4.2,  the  measured  Doppler  shift  must  be  adjusted  using  the 
estimated  receiver  velocity  component  in  the  direction  of  the 
slant  range,  and  the  resulting  lines  of  position  (from  successive 
Doppler  measurements)  must  be  subsequently  adjusted  using  the 
receiver  velocity  in  the  direction  orthogonal  to  the  slant  range. 
Any  error  in  the  estimates  of  these  receiver  velocity  components 
will  therefore  lead  to  a  position  fix  error.  Two  different 
effects  arise  from  the  velocity  error  parallel  to  the  slant  range 
(due  to  a  constant  error  and  a  differential  error)  which  are 
explained  in  detail  in  Sections  6.1  and  6.2  below.  Two  other 
effects  arise  from  the  velocity  error  orthogonal  to  the  slant 
range,  as  described  in  detail  in  Section  6.3.  Section  6.4  then 
summarizes  the  combined  effect  of  the  velocity  error  on  Transit 
position. 


6.1  Velocity  Error  Parallel  to  the  Slant  Range 


The  first  error  will  be  called  Arl,  and  is  due  to  the 
effect  of  a  constant  velocity  error  component  in  the  direction  of 
the  slant  range  p  from  the  satellite  to  the  receiver.  This 
component  of  the  velocity  error  directly  effects  the  Doppler 
measurement,  and  therefore  introduces  an  error  in  the  calculation 
of  the  cone  angle  y  (between  the  slant  range  and  the  satellite 
velocity  vector)  as  found  from  equations  (10)  and  (11). 

In  this  subsection  we  will  find  the  relationship  between 
this  velocity  error  and  the  resulting  position  error,  as  in 
equation  (2)  (assuming  zero  height  error).  To  do  this  we  first 
express  the  position  error  (AN,AE)  as  a  function  of  the  Doppler 
cone  angle  error  Ay,  and  then  we  express  the  Doppler  cone  angle 
error  as  a  function  of  the  velocity  error  (Avn,Ave). 


6.1.1  Magnitude  of  tu rl(Ay) 

In  expressing  the  position  error  Arl  as  a  function  of  the 
cone  angle  error  Ay,  we  first  find  its  magnitude,  then  its 
direction. 

Figure  6a  shows  how  a  constant  cone  angle  bias  Ay  (true  - 
estimate)  causes  the  position  estimate  to  shift  in  the  direction 
orthogonal  to  the  slant  range  vector  p  in  the  locally  horizontal 
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plane.  From  this  we  can  see  that  for  small  errors  Ay,  the 
position  error  vector 


Arl  ■  r  -  i 


(18) 


has  magnitude 


Arl  ■  |  Arl  |  *  p  sin|Ay| 


(19) 


6.1.2  Direction  of  Arl 

Figure  6a  shows  how  a  constant  cone  angle  bias  Ay  (true  - 
estimate)  causes  the  position  estimate  to  shift  in  the  direction 
orthogonal  to  the  slant  range  vector  p  (in  the  local  horizontal 
plane  at  the  receiver). 

Figure  6b  shows  the  situation  at  the  receiver  in  more 
detail.  Since  the  bearing  to  the  subpoint  is  the  horizontal 
projection  of  the  slant  range  vector  £  is  in  the  direction 
Thus  the  position  error  Arl,  which  is  orthogonal  to  the  direction 
i|>,  is  either  in  the  direction  (t|»n/2)  or  (^n/2). 

Thus  we  can  see  that  the  north  and  east  components  of 
this  position  error  are: 


'AIT|  Arl  cos('|%ll/2)  " 

.AeJ1  Arl  sin(*ten/2)  . 


(20) 


where  the  actual  sign  of  ±  is  yet  to  be  determined. 

At  the  point  of  closest  approach  p  is  also  orthogonal  to 
the  satellite  velocity  vector  V  (which  is  approximately  northward 
or  southward,  depending  on  the  sign  of  the  parameter  NS).  Again 
from  Figure  6a  we  can  see  that  the  error  Arl  (true  -  estimate) 
induced  by  a  positive  constant  cone  angle  bias  Ay,  is  roughly  in 
the  direction  -V.  Similarly  a  negative  Ay  would  cause  an  error 
Arl  in  the  general  direction  of  V. 
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North 


The  actual  direction  of  Arl,  as  seen  from  Figure  6b  is 
(<J»-ji/2)  rather  than  -V  and  (i(*ji/2)  rather  than  V.  This  is  the 
case  where  the  satellite  is  vest  of  the  receiver.  The  situation 
is  the  reverse  if  the  satellite  is  east  of  the  receiver:  the 
direction  of  Arl  is  (iJh-ji/2)  when  Ay  is  positive,  and  (^-n/2)  when 
Ay  is  negative.  This  east-vest  relationship  between  the  satellite 
and  the  receiver  is  contained  in  the  parameter  EV,  which  is 
provided  by  the  receiver  as  part  of  the  position  fix  data. 

Thus  the  direction  of  Arl  (the  sign  of  n/2  in  (20)),  is 
determined  entirely  by  i)/,  the  sign  of  Ay,  the  direction  of  V,  and 
EV.  Appendix  0  provides  the  expressions  needed  to  find  + 
(actually  only  sini|/  and  cos\|/  are  needed).  The  sign  of  Ay  will  be 
derived  in  the  next  section,  and  for  now  can  be  expressed  as 
Ay/|Ay|.  The  direction  of  V  for  this  purpose  is  contained  in  the 
variable  NS  (which  is  +1  if  the  satellite  is  moving  southward  and 
-1  if  moving  northward). 

A  change  in  sign  of  NS  corresponds  to  a  change  in 
direction  of  V  and  hence  a  change  in  direction  of  Arl.  A  change 
in  sign  of  Ay  has  the  same  effect  on  Arl.  Therefore  we  can  write 


U  -  *  ±  NS-EV-Ay  n/2  (21) 

WT 


where  now  the  ±  has  a  constant  value. 

Careful  consideration  of  the  example  shown  in  Figure  6b 
indicates  that  the  correct  value  is  negative.  In  the  case 
illustrated  by  Figure  6b  a  receiver  velocity  error  Av,  with  a 
component  towards  the  satellite,  has  caused  a  positive  cone  angle 
error  Ay  (the  Doppler  cone  angle  estimated  is  too  small,  as  seen 
from  equation  11).  This  positive  Ay  moves  the  position  fix,  f,  in 
the  direction  »|h-ti/2.  This  corresponds  to  a  position  error  vector, 
Arl,  in  the  opposite  direction,  i|*-n/ 2. 

Thus  the  sign  in  (20)  should  be  negative  in  the  situation 
illustrated  by  Figure  6b.  In  this  case  the  satellite  is  on  the 
south-to-north  half  of  its  orbit  (NS«-1)  and  the  receiver  is  east 
of  the  satellite  subpoint  (EV*-1).  The  cone  angle  error  Ay  here 
is  positive.  Thus  the  sign  in  (21)  must  be  negative,  so  that  (20) 
can  be  rewritten: 


'AN' 

f  Arl  cos(*  -  NS-EV-Ay  n/2) 

sAE. 

1 

Arl  sin(*  -  NS-EV-Ay  n/2) 

1 

l  1  Ay  | 

(22) 
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Careful  use  of  reduction  formulae  reduces  (22)  to 


n 

-  NS-EV-Ay 

'  Arl  sini|»  ) 

IaeJ 

1  Ay  | 

„  -Arl  cosi|>  J 

1 

(23) 


In  the  last  subsection  ve  found  the  magnitude  of  the  position 
error  Arl,  given  by  equation  (19).  Substituting  this  into  (23) 
yields: 


'AN' 

-  NS-EV‘Ay 

’  p  sin  |  Ay|  sin\|>  ' 

.AE„ 

i  T^rT 

„  -p  sin  |Ay|  cost|/  , 

’AN' 

'  p  sinAY  sin^  ' 

-  NS-EV 

.AE. 

.  -p  sinAY  cos^  , 

(24) 


(25) 


It  remains  now  to  find  dr  as  a  function  of  the  velocity  error. 


6.1.3  Cone  Error  as  a  Function  of  Velocity  Error 


Nov  to  find  the  relationship  betveen  the  velocity  error 
(Avn,  Ave)  and  the  cone  angle  error,  Ay,  ve  use  equation  (10), 
vhich  relates  the  range  rate  dp/dt  to  the  cone  angle.  Ve  can  see 
hov  their  errors  are  related  by  differentiating  (10)  vith  respect 
to  ys 


3p 

3  Y 


V  sinY 


(26) 


so  that  for  small  Ap 
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V  siny  Ay 


(27) 


Ap 


At 


Ap 


V  siny 


(28) 


This  shows  the  expected  change  in  y  due  to  satellite 
motion  along  a  straight  line,  as  a  function  of  the  expected 
change  in  dp/dt  (a  purely  geometric  relationship). 

Now  the  error  (true  -  estimate)  in  range  rate  Ap  caused 
by  a  velocity  error  Av,  is  the  projection  of  the  velocity  error 
Av  along  the  slant  range  vector,  as  seen  from  Figure  6b  (where 
the  angle  v  is  in  the  locally  horizontal  plane,  and  a  is  in  a 
locally  vertical  plane): 


Ap  *  Av  costi 


(29) 


where  Av  ■  | Av |  is  the  magnitude  of  the  velocity  error. 

The  sign  in  (29)  arises  from  careful  consideration  of 
Figure  6b,  where  it  can  be  seen  that  when  Av  (true  -  estimate)  is 
towards  the  satellite  (ie.  cosh  is  positive),  then  the  magnitude 
of  the  range  rate  estimate  |d p/d 1 1  will  be  too  large.  Since  dp/dt 
is  negative,  this  implies  that  the  range  rate  error  A*  p  (true  - 
estimate)  is  positive. 

For  example,  if  the  true  velocity  were  zero  and  the 
estimated  velocity  was  away  from  the  satellite  (so  that  Av  is 
towards  the  satellite  as  in  Figure  6b),  then  the  estimated  range 
rate  would  be  incorrectly  compensated  to  give  it  a  larger 
negative  value  (to  counter  the  effect  of  the  estimated  velocity, 
which  would  have  decreased  the  Doppler  shift  since  it  is  away 
from  the  satellite).  This  therefore  leads  to  an  estimated  range 
rate  that  is  too  small  (or  too  large  in  a  negative  sense).  In 
this  case  the  range  rate  error  (true  -  estimate)  would  be 
positive. 

Since  o,  h  and  (2n-v)  form  a  right  spherical  triangle, 
with  n  the  hypotenuse,  equation  (29)  can  be  written: 


Ap  »  Av  cos®  cosv  (30) 
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Now  from  (28)  and  (30),  the  Doppler  error  can  be  written  as: 


6y 


Av  coso  cosv 
V  siny 


(31) 


Since  we  have  assumed  that  the  velocity  error  is  much  smaller 
than  the  satellite  velocity  (  Av  «  V  )  equation  (31)  produces  a 
small  Ay  angle  (if  we  keep  in  mind  that  near  the  point  of  closest 
approach  siny  «  1).  Thus  we  have  sinAy  =  Ay  and  from  (25)  and 
(31): 


AN' 

*  NS-EV  p  Av  coso  cosv 

'  sint  ' 

.AE. 

-  Vsiny 

. -cost  , 

(32) 


At  the  PCA  y*Jl/2,  so  that  siny=l,  and  this  simplifies  slightly. 
The  effect  of  the  velocity  components  is  resolved  through  the 
angle  v.  From  Figure  6b  we  see  that  (since  Aa  here  is  a  negative 
angle) 


v  *  t  -  Aa 

(33) 

cosv  »  cos(t-Aa) 

(34) 

>  cosAot  cost  + 

sinAa  sint 

(35) 

a  (  Avn  cost  + 

Ave  sint  )/Av 

(36) 

Substituting  this  into  equation  (32)  now  allows  us  to 
relate  the  position  error  components  AN  and  AE  to  the  receiver's 
north  and  east  velocity  error  components  Avn  and  Ave  at  the  PCA 
(where  siny*l)s 
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(AN' 

, 

NS-EV  p  cos a 

'  (  COSl|> 

Avn  + 

sint|/  Ave 

)sin^  ' 

Iae. 

V 

-(  cos^ 

Avn  + 

sinij/  Ave 

)cost|/  . 

(37) 


'£N' 

as 

NS-EV  pcosa 

'  cos^sin^ 

sin2\|/  ' 

'Avn' 

.AE. 

1 

V 

.  -COS2i|» 

-cos^sin^  . 

.Ave. 

(38) 


In  the  appendices  it  will  be  shown  how  p,  t  and  V  can  be 
expressed  as  functions  of  the  known  quantities  X,L,a,  NS  and  EV. 
Appendix  A  gives  p(cr)  and  6(a).  Appendices  A,  B,  C  and  D  together 
give  4<a,  X,NS,EV).  Appendices  A  and  C  with  equation  (5)  give 
V(  a,  X,NS,EV) . 


6.1.4  Qualitative  Verification 


Since  the  Transit  satellites  are  in  polar  orbits,  we  can 
see  from  Figure  1  that  i|»  is  generally  close  to  ±90°,  in  which 
case  |sin^|  »  |cos^|  and  the  dominant  term  in  equation  (38)  is 
the  sin2'!/  term,  which  implies  that  Arl  is  largely  in  the 
north-south  direction,  due  to  the  east  velocity  error  (as  would 
be  expected  from  the  physics). 

A  simple  thought  experiment  can  verify  that  the  sign  of 
this  dominant  term  in  equation  (38)  is  correct.  If  we  consider 
the  situation  where  : 

NS  -  EV  -  1 

vn  «  Vn  ■  0 

ve  »  |  e  | 


then 


Ave  >  |  e  | 
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In  this  case  the  satellite  is  moving  southward,  and  is 
east  of  the  receiver.  Since  the  true  velocity  here  is  towards  the 
satellite  (ie.  to  the  east),  the  measured  Doppler  shift  will 
greater  than  would  be  received  by  a  stationary  receiver.  This 
causes  a  smaller  cone  angle,  moving  the  position  estimate  south 
of  the  true  position.  Therefore  AN  will  be  positive  in  this 
situation.  Thus  a  positive  Ave  produces  a  positive  AN  when 
NSsEV=l,  which  is  consistent  with  equation  (38). 

Another  way  to  arrive  at  the  same  conclusion  in  this 
situation  is  to  consider  that  the  Doppler  shift  will  be  zero 
after  the  true  point  of  closest  approach  (in  the  earth  fixed 
locally  level  frame),  at  which  time  the  satellite  will  be  south 
of  the  true  PCA.  This  again  implies  that  the  position  estimate 
will  be  south  of  the  true  position,  so  that  AN  will  be  positive. 
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6.2  Effect  of  Differential  Cone  Angle  Error 


From  equation  (31)  we  can  see  how  the  cone  angle  error  Ay 
changes  as  a  function  of  y  and  v.  Figure  6a  showed  how  a  constant 
cone  angle  error  Ay  causes  a  cross-range  error  Arl.  Figure  7 
illustrates  how  a  change  in  this  cone  angle  error  of  &y  causes  an 
error  in  range  Ap,  which  projected  on  the  ground  becomes  an 
along-range  position  error  Ar2.  In  this  section  we  will  find  the 
expression  for  this  position  error.  As  before,  we  first  find  the 
magnitude  Ar2  and  then  the  direction. 


6.2.1  Magnitude  of  Ar2 


In  Figure  7  we  let  yl  and  y2  (where  y2  a  n/2  is  at  the 
apparent  point  of  closest  approach)  be  the  cone  angle  estimates, 
taken  a  short  time  apart  (dt  seconds).  From  Figure  7  we  can  see 
that  the  differential  cone  angle  error  Sy  (the  change  in  the 
error  over  the  interval  dt)  causes  an  error  in  the  slant  range 
Ap.  From  the  right  angle  triangles  in  Figure  7,  we  have: 


tanyl  (39) 

and 

tan(  yl  -  8y)  =  (40) 

where  the  slant  range  error  Ap  is  related  to  the  horizontal  error 
magnitude  Ar2  by: 

Ap  =  Ar2  cosa  (41) 

Substituting  (39)  and  (41)  into  (40)  yields: 

tan( yl  -  5y)  =  tanyl  -  ^^0S<T  (42) 
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vhich  can  be  written 


Ar2  cos  a 
Vdt 


tanyl  -  tan(yl  -  8y) 


sin( &y) 

cosy!  cos(yl  -  8y) 


sinyl  f  sin(  Sr) _  ' 

cosyl  sinyl  cos(yl-Sy)  , 

Prom  (39)  this  becomes: 


Ar2  coso 
Vdt 


P  sin(  &y) _  ' 

Vdt 

.  sinyl  cos(yl-Sy)  , 


Therefore  the  range  error  magnitude  can  be  expressed 


_ p  sinSy _ 

Ar2 

sinyl  cos(yl  -  Sy)  coso 


To  eliminate  yl  from  (47),  we  use  equation 
gives  us: 


Vdt 

P 


cot(  yl) 


tan(n/2  -  yl  ) 


(43) 

(44) 

(45) 

(46) 

as: 

(47) 

(39),  which 

(48) 
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Since  the  cone  angle  y  is  taken  from  the  satellite 
velocity  vector  to  the  slant  range  vector,  it  is  always 

increasing  in  time.  Therefore  yl  <  y2  ■  n/2,  and  the  denominator 
of  equation  (46)  cannot  be  zero.  For  small  dt,  (in  the 

neighbourhood  of  the  PCA)  yl  s  n/2,  so  that  (48)  can  be  written 


Vdt 

P 


n/2  -  yl 


(49) 


Thus  the  cone  angle  estimate  at  a  short  time,  dt,  before  the 
point  of  closest  approach  (ie.  for  small  dt)  is: 


(50) 


Thus,  to  eliminate  yl  in  (47)  we  keep  the  first  order  terms  in 
dt,  which  are: 


and 


sinyl 


cos(yl-Sy)  «  cos(n/2  -  Vdt/p  -  Sy) 


sin(Vdt/p  +  Sy) 


«  Vdt/p  +  Sy 


Using  this  in  (47)  gives: 


(51) 

(52) 

(53) 

(54) 


(55) 
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Nov  ve  need  the  differential  expression  for  $y  (as  a  function  of 
dt).  To  obtain  this  ve  use  (31),  and  recognize  that  v  also 
changes.  (Since  ve  are  at  the  maximum  elevation  angle,  de/dt  =  0, 
and  a  can  be  considered  constant.  It  can  also  be  shovn  that  this 
is  the  point  of  closest  approach,  so  that  p  is  at  a  minimum  and 
dp/dt  >0.  Thus  p  can  also  be  considered  constant.) 


5y  ■  Ay2  -  AyI 


(56) 


Av  COSIT 
V 


'  cosv2 
.  sinY2 


cosvl  ' 
sinyl  . 


(57) 


vhere  Av  is  the  velocity  error  vhich  induced  these  cone  angle 


errors.  Using  (50)  for  yl»  and  since  y2  -  n/2,  this  becomes: 


5y  2  -Avcoscr 


cosvl 


sin( n/2-Vdt/p) 


-  cosv2 


(58) 


-Av  cosq  (  cosvl 
V  lcos(Vdt/p) 


cosv2 


(59) 


From  Figure  8  ve  vill  see  that  the  angle  d£  betveen  slant  ranges, 
is  simply  related  to  the  angle  in  the  horizontal  plane: 


dv  a  v2  -  vl 


(60) 


vhich  is  the  difference  in  bearing  from  the  velocity  error  to  the 
horizontal  projection  of  the  slant  ranges  (as  can  be  seen  from 
Figure  6b).  Since  these  bearings  are  defined  by  the  geodesics 
(vhich  are  great  circles  if  a  spherical  earth  model  is  used)  from 
the  receiver  to  the  subpoints,  and  these  great  circles  are  the 
intersection  of  vertical  planes  vith  the  earth's  surface,  ve 
therefore  must  project  dovn  from  the  slant  ranges  using  vertical 
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Figure  8.  Vertical  Projection  From  Slant  Range  to  Horizontal. 


-34- 


planes,  as  shovn  in  Figure  8. 


Consider  the  small  right  spherical  triangle  MEF  in 
Figure  8.  A  trigonometric  identity  gives  the  relationship 


sin(n/2-a)  =  tan|dE|*cotA 

(61) 

vhich 

can  be  rearranged  as 

tan|d£|  =  coso’tanA 

(62) 

Nov, 

same 

by  considering  the  large  right  spherical 
trigonometric  identity  yields: 

triangle  AABD 

sin(n/2)  *  tan|dv|‘COtA 

(63) 

• 

•  « 

tan|dv|  =  tanA 

(64) 

From  (62)  and  (64)  we  therefore  have 

tan|dv|  *  tan  |d C I  (65) 

cos  a 

Therefore,  for  small  dv  ve  have  (and  since  0  <  c  <  n/2  ) 


|dv|  «  Jd£l 
cos  a 


(66) 


From  Figure  6b  ve  can  see  that  v  is  increasing  in  time  (dv  is 
positive)  if  the  satellite  is  to  the  east  moving  south  or  to  the 
vest  moving  north,  othervise  dv  is  negative.  In  other  vords,  the 
sign  of  dv  is  NS’EW.  Thus  equation  (66)  can  be  revritten: 


(67) 


dv  -  NS-BW-  ld£| 
cos  a 

From  Figure  7  ve  can  see  that 

|d£|  -  n/2  -  yl  (68) 

From  (50)  this  gives,  to  first  order  in  dt: 

|d£|  «  Vdt/p  (69) 


Thus,  from  (60),  (66)  and  (69)  ve  have: 


v2  *  vl 

+  NS-EV-  |d  ^ | 

(70) 

COSO 

a  Vl 

+  NS-EV-Vdt 

(71) 

pcos  a 

Therefore 


cosv2  a  cos /vl  +  NS’EV^Vdt  'j  (72) 

V  pcoso  ) 


»  cosvl  cos  /  Vdt  ) 

-  NS-EV-sinvl  sin  /  Vdt  ) 

(73) 

(pcos  a) 

(pcos  a) 

Nov  the  term  Vdt/p  is  very  small,  since  the  slant  range  p  is 
greater  than  the  orbital  height  (1,075,000  metres)  and  the 
satellite  velocity  V  is  less  than  7,310  m/s  (from  (7)).  Thus  for 
small  dt  ve  can  approximate,  to  first  order  in  dt: 


cosv2  «  cosvl  -  NS*EV*sinvl*_Vdt_  (74) 

pcos  o 
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This  can  now  be  used  in  (59)  to  express  8y  in  terms  of  dt, 
keeping  only  the  first  order  terms  in  dt,  (also,  since  ve  have 
eliminated  the  m2  terms,  ve  can  just  use  v  for  vl): 


8y  *  -Av  cosff  (  cosv  +  NS*EV*sinv  Vdt  (75) 

V  l  1  '  cosv  J 


(76) 


This  can  now  be  used  in  (55)  to  express  the  position  error 
independently  of  dt: 


Ar2  a  p  ( -NS •  EV  Av  sinv  dt/p) _  (77) 

[(-NS’EV  Av  sinv  dt/p)  +  Vdt/pjcoso 


-NS-EW  p  Av  sinv _  (78) 

-NS*EV  Av  coso  sinv  +  Vcoso 


Nov  for  velocity  errors  Av  much  smaller  than  the  satellite 
velocity  V,  this  can  be  vritten: 


(79) 


This  can  be  expressed  in  terms  of  the  components  of  the  velocity 
error  vector  Av  by  using  the  fact  that  v  -  *  -  Aa  (see  equation 
(33)),  vhere  Aa  is  the  bearing  of  Av. 
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sinv 


sin('^-Aa) 


(80) 


-  sin^  cosAa  -  cos\p  sinAa  (81) 

Av  sinv  -  simj<Av  cosAa)  -  cosiKAv  sinAa)  (82) 
■  siniji  Avn  -  costy  Ave  (83) 

Thus  from  (79)  we  have: 

Ar2  ■  -NS- EV  p  (sim|*  Avn  -  cos^  Ave)  (84) 

Vcosff 


6.2.2  Direction  of  Ar2 


Nov  ve  must  express  the  components  of  the  position  error 
(AN, AE)  in  terms  of  the  position  error  magnitude  Ar2.  From  Figure 
9  it  can  be  seen  that  a  positive  differential  cone  angle  error  Sy 
(defined  by  (56))  causes  the  position  estimate  to  move  in  the 
direction  avay  from  the  subpoint  (at  bearing  \Hn).  Thus  a 
positive  &y  creates  a  position  error  vector  Ar2  ( true-estimate) 
with  bearing  From  (76)  and  (79)  ve  can  see  that  &y  has  the 
same  sign  as  Ar2  (both  are  determined  by  the  sign  of  sinv). 
Therefore  Ar2  has  bearing  t|*Jt  if  Sr  (given  by  (76))  is  positive, 
and  has  bearing  ^  if  8y  is  negative.  Since 


cos (^h- it)  -  -cos^ 
sin(i|H-n)  »  -sinij/ 
ve  can  therefore  vrite: 


(AN] 

'  cosy  ' 

* 

Ar2 

IaeJ 

>  siniji  . 

(85) 


(86) 
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Figure  9.  Direction  of  A  r 2. 


<«- 


Finally  we  can  use  (84)  to  substitute  for  Ar2  in  (86)  to  obtains 


NS-EW  p 
Vcoso 


(sintf;  Avn 
(sin\|/  Avn 


cos\|/  Ave)  cost); 
cost);  Ave)  sint); 


J'AN'l 

'  -sint)cost); 

COS  2  tj;  > 

I'Avn'j 

a 

NS-EV  p 

UJ2 

Vcoser 

k  -sin*t); 

costjsint);  „ 

UveJ 
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6.3  Velocity  Error  Orthogonal  to  the  Slant  Range 


The  remaining  velocity  induced  position  errors  are  due  to 
the  effect  of  the  velocity  error  component  orthogonal  to  the 
slant  range.  In  Section  6.3.1  ve  will  see  how  this  affects  the 
individual  cone  angle  calculations,  and  in  Section  6.3.2  we  will 
see  how  it  causes  an  additional  error  when  cones  are  translated 
to  remove  the  effect  of  receiver  motion  between  successive  cone 
angle  measurements  (prior  to  intersecting  the  cones). 


6.3.1  Cone  Angle  Error 


A  velocity  error  Av  which  is  orthogonal  to  the  slant 
range  will  be  parallel  to  the  satellite  velocity  vector  V  (at  the 
point  of  closest  approach).  This  will  not  affect  the  cone  angle 
at  the  point  of  closest  approach  (as  seen  in  Section  6.1  above), 
but  it  will  effect  the  cone  angle  measurement  on  either  side  of 
the  PCA,  where  there  will  be  a  component  of  Av  along  the  slant 
range.  Hore  importantly,  this  component  of  Av  will  be  positive  on 
one  side  of  the  PCA  and  negative  on  the  other,  which  will  cause 
the  cone  angle  to  increase  on  one  side  of  the  PCA  and  decrease  on 
the  other.  This  will  move  the  intersection  of  the  cones  towards 
or  away  from  the  satellite,  introducing  a  position  error.  The 
situation  is  illustrated  in  Figure  10. 


6. 3. 1.1  Magnitude  of  Ar3 


Figure  10  illustrates  the  situation  where  the  satellite 
is  at  SI,  a  short  time  dt  before  the  point  of  closest  approach 
(PCA),  and  therefore  is  at  a  distance  of  Vdt  from  the  PCA.  We  let 
D  be  the  slant  range  at  SI,  and  r  be  the  cone  angle.  Then  the 
cone  angle  error  Ay  at  SI  causes  a  slant  range  error  Ap,  as 
shown.  Simple  trigonometry  gives  us: 


Ap  cos(  r+Ay)  =>  D  sin  Ay 


(89) 


and 


D 


p/siny 


(90) 


so  that 
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Ap  =  p  sinAy  (91) 

sin-ycos(r+Ar) 


Nov  equation  (31)  of  Section  6.1.3  gives  us  the  expression  for 
the  cone  angle  error  Ay  caused  by  a  velocity  error  Av: 


Ay  =  Av  cosg  cos£ 
V  sinv 


(92) 


where  a  is  the  satellite  elevation  angle  and  y  is  the  cone  angle 
at  Si.  Here  ve  use  £  for  the  angle  betveen  the  velocity  error  Av 
and  the  bearing  to  the  subpoint  for  SI,  so  that  ve  can  reserve  v 
(as  shown  in  Figure  6b)  for  the  angle  betveen  Av  and  the  bearing 
line  to  the  subpoint  at  PCA. 

From  equations  (67)  and  (69)  ve  saw  that  the  bearing 
angle  £  a  v  -  dv,  a  time  dt  before  the  PCA,  is  related  to  the 
angle  v  at  the  PCA  by 


£  «  v  -  NS-EV  Vdt 
pcos  a 


(93) 


so  that 


cos£  a  cosv  cos  f  Vdt 

+  NS’EV-sinv  sin  | 

'  Vdt  \ 

(94) 

Ipcos  a) 

1 

vpcos  a) 

=  COSV  +  NS'EV 

•sinv*  f  Vdt  1 
(pcos  oj 

(95) 

Using  this  in  (92)  gives: 

Ay  -  Av  cos  «x  (cosv 
VsinY  l 

+  NS-EV-sinv  Vdtl 
pCOSff  J 

(96) 

*  Av  cos  a  cosv 

NS'EV’  Av  sinv  dt 

(97) 

VsinY 

p  sinY 

Nov  the  first  term  is  from 

the  component  of  Av  parallel  to 
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slant  range  (Avcosv),  and  is  symmetric  in  the  cone  angle  v.  This 
gives  rise  to  a  position  error  orthogonal  to  the  slant  range, 
(and  in  fact  is  the  error  Ay  described  in  equation  (31),  which 
led  to  Arl),  which  does  not  contribute  to  the  error  described  by 
Figure  10.  The  second  term  however,  is  antisymmetric  in  the  cone 
angle  v  (  since  sin(-v)  -  -sinv  ),  causing  an  along-range  error 
which  we  call  Ar3. 

Using  this  second  term  of  equation  (97)  in  (91),  (and 
since  Ay  is  a  small  angle:  sinAysAy  and  cos(  y+Ay)=cosy  ),  we 
obtain 


Ap  =  NS^EW*  Av  sinv  dt  (98) 

sin2y  cosy 

From  equation  (50)  we  saw  that  near  the  PCA 

y  «  n/2  -  Vdt/p  (99) 

so  that,  since  Vdt/p  «1 


siny  =  1 
cosy  a  Vdt/p 


(100) 


Using  this  in  (98)  gives 


Ap  =  NS'EV-p  Av  sinv  (101) 

V 


Projecting  this  into  the  horizontal  plane,  as  illustrated  in 
Figure  7  for  Ar2,  we  obtain  the  magnitude  of  Ar3: 


Ar3  =  NS*EW*p  Av  sinv 
V  cosc 


(102) 


Note  that  this  position  error  is  independent  of  the  time  dt 
before  (or  after)  the  point  of  closest  approach,  so  that  the  cone 
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error  described  by  (97)  does  lead  to  a  constant  position  error. 

Furthermore,  this  error  Ar3  can  be  seen  to  have  exactly 
the  same  magnitude  as  Ar2  (by  comparing  equation  (102)  to 
equation  (79)).  Ve  could  therefore  similarly  express  Ar3  in  terms 
of  +  rather  than  v  by  folloving  the  same  derivation,  given  by 
equations  (79)  to  (84). 


6. 3. 1.2  Direction  of  Ar3 


As  shown  in  Figure  10,  the  direction  of  Ar3  is  towards 
the  subpoint  (bearing  <10  if  the  velocity  error  Av  has  a  positive 
component  in  the  satellite  direction  of  travel  V.  Conversely  Ar3 
is  away  from  the  subpoint  if  Av  is  towards  -V.  Thus  Ar3  can  be 
expressed  as  in  equation  (88),  with  a  possible  sign  change.  Using 
Figure  10  to  test  the  sign,  we  see  that  with  NS  *  EV  =  -1,  a 
positive  Avn  should  cause  a  negative  AE,  which  equation  (88) 
does.  Therefore  Ar3  is  identical  to  Ar2,  and  is  therefore  given 
bys 


'AN' 

* 

NS-EW  p 

'  -sinij/cosijf 

cos2V*  ' 

'Avn' 

.AE. 

'X 

Vcoso 

„  -sin*i|> 

cos\|einij/  . 

.Ave, 

(103) 
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6.3.2  Translation  Error 


Since  the  receiver  moves  between  Doppler  measurements, 
the  Doppler  cones  must  be  translated  by  the  displacement  vector 
(v  dt)  in  the  direction  of  the  receiver  velocity  to  form  these 
intersections.  Figure  11a,  shoving  the  "slant  range  plane", 
illustrates  how  an  error  in  this  velocity  Av  over  an  interval  At, 
moves  the  cone  by  an  excess  amount  (Av  dt)  which  moves  the  cone 
intersection  in  the  direction  of  the  slant  range  by  an  excess 
amount  Ap.  The  resulting  position  error  Ar4  is  the  projection  of 
this  slant  range  error  Ap  into  the  local  horizontal  plane 
(perpendicular  to  the  slant  range,  as  shown  in  Figure  7).  The 
position  r  and  its  estimate  are  more  properly  shown  in  the 
horizontal  plane  as  in  Figure  lib,  which  shows  the  quantities  of 
interest,  assuming  zero  true  velocity  for  simplicity. 

In  the  receiver's  horizontal  plane,  as  shown  in  Figure 
lib,  the  cones  (on  intersection  with  the  horizontal  plane)  reduce 
to  lines  of  position  (LOPs).  In  this  plane,  in  the  neighbourhood 
of  the  receiver,  these  LOPs  are  just  the  projections  of  the  slant 
range  vectors  into  the  horizontal  plane,  as  shown  in  Figure  12. 
This  projection  is  in  the  direction  of  the  line  of  intersection 
of  the  two  cones  (orthogonal  to  both  the  slant  range  p  and  the 
satellite  velocity  V,  as  seen  in  Figure  2). 


6. 3. 2.1  Magnitude  of  Ar4 


From  Figure  lib  we  can  see  that  the  magnitude  T  of  the 
component  of  the  velocity  error  (Av  dt)  orthogonal  to  the  LOP  is 
found  by  projecting  with  the  included  angle  \>  (where  v  is  defined 
as  shown  in  Figure  6b).  Thus  the  effective  LOP  displacement  T 
(orthogonal  to  the  LOP)  due  to  this  velocity  error  has  magnitude: 


f  =  Av  dt  |sinv| 


(104) 


Now  the  resulting  position  error  Ar4  can  be  seen  from 
Figure  lib  to  be  parallel  to  the  slant  range  and  of  magnitude  Ar4 
related  to  the  displacement  error  by: 


Ar4  |sin(d£)  |  =  T 


(105) 


where  dC  is  the  change  in  the  angle  between  lines  of  position 
(intersection  of  Doppler  cones  with  the  horizontal  plane)  over 
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the  time  interval  dt  (assuming  a  constant  velocity  error  Av  over 
this  interval).  From  (104)  and  (105)  ve  have: 


Ar4  |sin(dC)  |  ■  Av  dt  |sinv|  (106) 


For  small  intervals  dt,  clearly  dC  is  also  very  small,  so  that  ve 
can  use  sin(dC)  =  dC,  and  vrite  (106)  as 


Ar4 


Av  |sinv| 

(TdtRdt) 


(107) 


which  relates  the  velocity  error  magnitude  Av  to  the  magnitude  of 
the  position  error  Ar4. 

Consider  the  spherical  triangle  AGBD  in  Figure  12.  This 
spherical  triangle  is  formed  by  approximating  the  Doppler  cones 
at  the  receiver,  near  the  point  of  closest  approach,  as  planes 
(tangent  to  the  cone  along  the  slant  range).  This  is  valid  since 
the  cone  at  the  PCA  does  reduce  to  a  plane,  and  in  any  case  the 
curvature  will  be  very  small  compared  to  the  position  error  Ar4. 
Note  that  since  the  circle  of  intersection  of  the  Doppler  cones, 
C,  is  orthogonal  to  the  slant  range  p: 


U  -  d£ 


(108) 


Note  also  that  the  smaller  spherical  triangle  AGEF  is  a 
hemisphere  slice,  so  that  its  angles  E  and  F  are  n/2  and  the 
sides  GB  and  GF  subtend  n/2  at  the  centre  r.  Therefore  the  lav  of 
cosines  applied  to  AGEF  gives  us 


cos(dQ  -  cos2(n/2+a)  +  sin2(n/2+o)  cosy  (109) 

-  sin2a  +  cos2o  cos(dS)  (110) 

Thus  for  small  dC  (using  cose  «  l-e2/2): 

1  -  (dOJ/2  -  sin2 a  +  cos2a  (  1  -  (d£)2/2  )  (111) 
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(dC)2 


cos2o  (d£)2 


(112) 


and  so 


|dt|  -  |d£|coso 


(113) 


Using  (69)  for  |d£|  this  becomes 


|d  C  |  «  Vdt  coso 

P 


(114) 


This  can  now  be  used  in  (107)  to  express  the  relationship  between 
the  position  error  magnitude  Ar4  and  the  velocity  error  Av  in 
terms  of  "knovn"  quantities. 


Ar4  «  Av  |sinvl  p 
Vcoso 


(115) 


6. 3. 2. 2  Direction  of  Ar4 


Nov  we  must  express  the  components  of  the  position  error 
( AN, AE)  in  terms  of  the  position  error  magnitude  Ar4.  From 
careful  consideration  of  Figures  11a  and  lib,  it  can  be  seen  that 
in  general  Ar4  is  in  the  direction  towards  the  subpoint  (bearing 
i|>)  if  the  projection  of  the  velocity  error  Av  in  the  satellite 
direction  of  travel  V  is  positive,  and  conversely  that  Ar4  is 
away  from  the  subpoint  (bearing  ij»n)  if  this  projection  is 
negative.  Since 


cos(t*-n)  ■  -cos<|# 
sin('HH)  ■  -sin<|» 


(116) 


ve  can  therefore  vrite: 
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rAN' 

'  COSl|f  ' 

-  ±Ar4 

.AE, 

A 

.  sin*  . 

(117) 


or,  using  equation  (115) 


'AN' 

_ 

±Av  sinv  p 

'  cosy  " 

4 

Vcoscr 

.  sini|/  , 

(118) 


where  the  ±  has  absorbed  the  absolute  value  sign  from  (115),  and 
hence  the  ±  of  equation  (118)  may  not  be  the  same  as  in  equation 
(117). 

In  the  particular  situation  shown  in  Figure  11a  (  NS— 1, 
EV— 1,  sin\K0  ),  the  velocity  error  Av  has  a  positive  component 
in  the  satellite  direction  of  travel  V  (V  and  Av  are  both  to  the 
north),  so  that  the  resulting  position  error  Ar4  has  bearing  ^ 
(towards  the  satellite  subpoint).  Therefore  the  sign  of  the  first 
term  in  equation  (118)  should  be  positive,  to  produce  a  vector  in 
the  direction  of  i|>.  Since  sinv  in  this  case  is  negative,  the  ±  in 
equation  (118)  should  also  be  negative. 

To  generalize  this:  there  are  two  factors  which  determine 
whether  or  not  Av  has  a  component  in  the  direction  of  V.  These 
are: 


-the  direction  of  V,  as  determined  by  NS 
-the  direction  of  Av 


The  second  factor  here  can  in  turn  be  determined  by  two  factors: 

-the  angle  v  between  Av  and  the  subpoint  vector, 

-the  direction  to  subpoint  vector,  as  determined  by  EV 


Ve  can  see  therefore  that  the  sign  in  equation  (118) 
should  change  with  NS,  since  the  direction  of  V  does.  Ve  can  also 
see  that  the  sign  in  equation  (118)  should  change  with  the 
relative  direction  of  Av,  which  can  be  accounted  for  by  the  sign 
of  sinv  and  by  EV.  Ve  can  therefore  express  (118)  as: 
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'AN' 

=  -NS-BV  p  Av  sinv 

'  cost|»  " 

.AE> 

.  Vcoso 

.  sin $  t 

(119) 


Finally  ve  can  use  (83)  to  substitute  for  Av  sinv  in  (119)  to 
obtain: 


AN' 

■  -NS'EW  p  (simp  Avn  -  cos^  Ave) 

'  cos^  " 

.  Vcosct 

.  sin*  „ 

(120) 


or 


AN' 

NS-EW  p 

'  -sin^cos^ 

cos2^  ' 

Avn' 

M, 

4 

Vc  os  a 

>  -sin2^ 

cos^sin^  , 

>Avey 

(121) 


Comparing  this  to  (88)  and  (103)  we  can  see  that  this  error  is  in 
fact  identical  in  magnitude  and  direction  to  Ar2  and  Ar3.  If 
there  is  a  simple  reason  for  this  it  is  not  immediately  apparent. 
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6.4  Coabined  Effect  of  Velocity  Error 


Equations  (38),  (88),  (103)  and  (121)  can  be  combined  to 
express  the  total  effect  of  the  error  in  the  velocity  provided  to 
a  Transit  receiver  on  the  Transit  position  fix  error: 

AN  *  (  l  coso  -  3/coso  ]cos\)sin\|»  dvn 

+  [  cososin2^  +  (3/coso)cos2t|>  Jdve  1  (122) 


dE  «  [  —cos ocos 2 ^  -  (3/coso)sin2^  Jdvn 


+ 


-coso  +  3/coso  Icosijsinij/  five 


) 


(123) 


or 


dN  *  NS-EWpf  -(3-cos2 o)cost|eini|>  dvn 
Vcosov 


+ 


cos2osin2i(> 


+ 


3cos 2 


(124) 


dE 


NS-EVp /  [ -cos 2  ocos 2 1|»  -  3sin2\|>  ]dvn 
Vcosol 

+  (3-cos2 o)cos\jein\|»  dve 


(125) 


Now  if  we  let 

A  ■  cos2o 
B  ■  3  -  cos2o 

then 


(126) 

(127) 
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(128) 


NS-Effp 

Vcosff 


-Bcos^slntfr  Asin2^f3cos2^r  ' 
-Acos2^-3sin2t|»  Bcos^sinf 


Together  with  equation  (17)  for  the  height  sensitivity 
components,  this  gives  us  all  components  of  the  Transit  position 
measurement  sensitivity  matrix  H,  of  equation  (1),  as  needed  by  a 
Kalman  filter  to  properly  process  a  Transit  fix.  These  components 
are  given  here  as  functions  of  p,  a,  \J/,  NS,  EV  and  V.  The 
receiver  generally  provides  a,  NS  and  EV  (as  well  as  X  and  L). 
The  remaining  variables  p,  t|/  and  V  can  be  expressed  as  functions 
of  these  receiver  supplied  quantities. 

As  mentioned  earlier,  it  is  noteworthy  that  in  most  cases 
the  bearing  to  the  subpoint  ij/  is  approximately  ±n/ 2,  so  that  cos\|» 
is  of  small  magnitude  relative  to  sintfr.  Note  also  that  at  high 
satellite  elevation  angles  (o*n/2)  the  1/coso  factor  becomes  very 
large. 


7.0  COMBINED  EFFECT  OF  VELOCITY  AND  HEIGHT  ERRORS 

Equations  (17)  and  (128)  can  be  combined  to  express  the 
total  effect  of  the  error  in  the  velocity  and  height  provided  to 
a  Transit  receiver  on  the  Transit  position  fix  error: 


H 


f-cost>»tano  -NS'EV  pBcos'lsiniJ; 
Vcoso 


NS-EVp  (Asin2t|»3cos2»)  'j 
Vcoso  j 


-sinij/tano  -NS-EW  p(Acos2^f3sin2\<;)  NSj_EWpBcos\|sin>j/ 
^  Vcoso  Vcoso 


(129) 


where 


A  ■  cos  2o 


(130) 


B  a  3  -  cos2o 


(131) 


so  that 


-C-BcoS'jsinV'  C(Asin2\|M-3cos2ij/) 

fW 

-C(Acos2*(>f3sin2i|0  C-Bcos^sin^ 

(132) 

where 

C  a  NS-EVp  (133) 

Vcoso 


This  provides  all  of  the  components  of  the  Transit 
position  measurement  sensitivity  matrix  H,  of  equation  (1),  as 
needed  by  a  Kalman  filter  to  properly  process  a  Transit  fix. 
Notice  from  (133)  that  C  diverges  as  the  elevation  angle  o 
approaches  90°  (while  A  goes  to  0  and  B  goes  to  3),  reflecting 
the  fact  that  high  elevation  angle  satellite  passes  cannot 
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provide  a  useful  fix  (due  to  Ar2,  &r3  and  Ar4). 

The  H  components  are  given  here  as  functions  of  p,  a, 

NS,  EV  and  V.  The  receiver  generally  provides  a,  NS  and  EV  (as 
veil  as  X).  The  remaining  variables  p,  +  and  V  can  be  expressed 
as  functions  of  these  receiver  supplied  quantities. 

Section  8  belov  provides  an  algorithm  by  which  these  H 
components  can  be  evaluated. 
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8.0  ALGORITHM 


In  this  section  ve  assemble  all  of  the  relevant  equations 
needed  to  calculate  the  H  matrix  from  the  four  input  quantities 
that  are  assumed  to  be  known: 

1/  receiver  latitude  X 
2/  satellite  maximum  elevation  angle  a 
3/  satellite  north- to-south  direction  of  travel  NS 
4/  satellite  east-to-west  direction  of  travel  EW 

Since  these  are  needed  after  the  fix  is  made,  the  fix  latitude 
can  be  used  for  X,  and  modern  receivers  usually  do  supply  o,  NS 
and  EV. 


The  following  set  of  equations,  gathered  from  throughout 
this  paper  and  its  appendices,  have  all  been  derived  herein  and 
are  sufficient  for  evaluating  the  components  of  H  as  functions  of 
(X,  <r,NS,EV).  These  equations  can  be  solved  in  the  order  presented 
here.  It  is  also  important  to  apply  the  consistency  check  (A15), 
which  limits  6  at  high  latitudes  X.  If  this  limit  is  violated 
(because  of  an  elevation  angle  a  that  is  too  low  for  the  given 
latitude  X)  then  the  square  root  in  (C15)  may  not  have  a  real 
value. 


(A12) 


(A15) 


(A5) 


K.  ■  (R+H)w/Vg  =  .0744716  (Cll) 


-58- 


-cos»(rtan<T  -NS’EVpBcosygin^  NS^EWp_(AsinJ\K3cos2i</)' 

H  »  Vcosa  Vcoso 

(129) 

-sln\)A:an0  -NSJSWp(Acos2 \j*3sin2 NS,EWpBcos4sin\t> 
w  Vcosa  Vcos a 


A  ■  cos 2  a  (130) 

B  ■  3  -  cos 2  a  (131) 
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The  tree  structure  shown  in  Figure  13  illustrates 
graphically  how  the  these  equations  can  be  used  to  find  the 
components  of  H.  Each  node  in  Figure  13  shows  the  equation  number 
followed  by  the  variable  which  that  equation  expresses  as  a 
function  of  known  quantities  and  the  unknown  variables  attached 
to  the  node  from  below.  Thus  for  example  equation  B13  expresses  $ 
in  terms  of  known  quantities  and  the  unknowns  X  and  EV.  To 
evaluate  any  variable  in  this  figure  one  must  therefore  solve  for 
all  variables  below  its  node. 


Figure  13.  Appendix  Use  Flow  Diagram. 


The  following  is  a  pseudo  code  algorithm  which  uses  these 
equations  to  evaluate  the  components  of  the  measurement  matrix  H, 
given  the  Transit  satellite  pass  information  (arNS,EW)  and  a 
receiver  latitude  estimate  X. 
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input  variables: 
output  variables: 
constants: 


Equation 


A,NS,EV, a 
H<2,3) 

R  -  6370000. 

H  =  1075000. 

Vs  -  7290 

K  a  .0744716  «(R+H)/Vs 

K1  -  .855608 
n  -  3.14159265 
intermediate  variables:  p,V,cscphi 

csig, c2sig,s2sig,  clat,slat, 
cthe,sthe,s2the,  cpsi,spsi 
dummy  variables:  A,B,C,D 


if(  a  <  .1  or  a  >  1.4  )  return  error  condition 

csig  *  cos  (a) 
c2sig  *  csig*csig 
s2sig  »  1  -  c2sig 

cthe  ■  Kl*c2sig  +  sqrt(  s2sig*(l  -  Kl*Kl*c2sig)  ) 

if(  cthe  <  abs(sin(A))  )  return  error  condition 

p  -  (  (R+B)*cthe  -  R  )/sin(<x) 

s2the  =  1  -  cthe*cthe 
sthe  =  sqrt(s2the) 

D  *  sqrt(  cthe*cthe  -  4*K*(NS*EV*sthe*sin(X)  -  K*s2the)  ) 
slat  -  (cthe  -  D)/(2*NS*EV*K*sthe) 
clat  a  sqrt(l  -  slat*slat) 

D  a  K*clat 

cscphi  a  -EV*sqrt(l  +  D*D) 

spsi  a  -clat/(cos(X)*cscphi) 

cpsi  a  (slat  -  cthe*sin( A))/(sthe*cos( A)) 

V  a  Vs*sqrt(l  +  K*K*clat*clat) 

D  a  sin(a)/csig 
HH(1, 1)  -  -cpsi*D 
HH(2,1)  a  -spsi*D 

A  a  c2sig 

B  a  3  -  c2sig 

C  -  NS*EV*p/(V*csig) 

HH(1,2)  ■  -C*B*cpsi*spsi 

HH(2,2)  -  -C*(A*cpsi*cpsi  +  3*spsi*spsi) 

HH(1,3)  -  C*(A*spsi*spsi  +  3*cpsi*cpsi) 

HH(2,3)  -  C*B*cpsi*spsi 
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(Cll) 

R/(R+H) 


(A12) 

(A15) 

(A5) 

(C15) 

(B13) 

(D3) 

(D6) 

(5) 

(17) 

(130) 

(131) 
(133) 


9.0  NUMERICAL  RESULTS 


From  Section  8.0  above,  we  see  that  the  Transit  error 
sensitivity  model,  as  embodied  in  the  H  matrix  of  (129),  is  a 
function  of  four  variables:  ( cr,NS,EV,  X) .  This  model  can  therefore 
be  illustrated  by  shoving  the  numerical  values  of  the  H  matrix  as 
functions  of  the  satellite  elevation  angle  or,  for  each  of  the 
four  possible  directions  of  travel  (NS=±1,EW=±1) ,  at  a  fixed 
receiver  latitude  X.  Section  8.0  above  provides  the  explicit 
equations  necessary  to  do  this.  In  Section  9.2  below  the  six 
elements  of  H  are  plotted,  along  with  three  intermediate 
geometric  quantities  of  interest,  for  each  of  the  four  directions 
of  travel,  at  a  receiver  latitude  of  45°.  The  latitude  dependence 
is  illustrated  in  Section  9.3,  by  considering  the  receiver 
latitude  to  be  a  parameter  defining  a  family  of  curves  for  each 
of  the  curves  shown  in  Section  9.2. 

In  Section  9.1  it  is  shown  that,  except  for  a  change  in 
sign,  there  are  actually  only  four  different  height  error 
sensitivity  functions  (out  of  a  potential  8,  from  2  H  elements 
times  4  directions  of  travel),  and  there  are  only  six  different 
velocity  error  sensitivity  functions  (out  of  the  potential  16, 
from  four  H  elements  times  4  directions  of  travel).  From  the 
plots  of  Section  9.3  it  can  be  seen  that  in  fact  one  of  these 
four  different  height  sensitivity  functions  (clockwise  H(2,l))  is 
actually  very  close  to  an  other  (counter-clockwise  H(2,l)),  and 
two  of  the  velocity  sensitivity  functions  (clockwise  H(2,2)  and 
H(l,3))  are  very  close  to  two  others,  (counter-clockwise  H(2,2) 
and  H(l,3))  so  that  there  are  in  effect  only  three  "basic"  height 
sensitivity  functions  and  four  "basic"  velocity  error  sensitivity 
functions  (except  for  a  change  in  sign). 


9.1  Direction  of  Travel  Dependence 


By  inspecting  the  equations  of  Section  8.0  above,  which 
express  the  H  components  explicitly  in  terms  of  the  variables 
( <t,NS,EV, Xp) ,  it  can  be  seen  that  the  direction  of  travel 
variables,  NE  and  EW,  appear  only  in  equations  (C15),  (B13)  and 
(107).  In  two  of  these  equations  they  appear  only  as  the  product 
NE-EW,  which  has  only  two  possible  values.  These  two  cases  can 
conveniently  be  labelled  the  "clockwise"  case,  for  NS'EW=1  (i.e. 
NS»EW*1  and  NS»EV=-1)  and  the  "counterclockwise"  case,  for 
NS-EV—1  (ie.  NS*1,EW=-1  and  NS=-1,EV=1),  describing  the 
satellite  direction  relative  to  the  receiver  looking  down  from 
above.  In  the  remaining  equation  (B13),  only  EW  appears,  and  its 
effect  is  simply  to  change  the  sign  of  ♦,  the  effect  of  which  can 
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be  easily  folloved  through  to  see  that  it  simply  changes  the  sign 
of  ^  (via  (D3)),  which  in  turn  changes  the  sign  of  H(2,l),  H(l,2) 
and  H(2,3). 

Therefore  ve  can  conclude  that  all  terms  for  NSsEV«l  will 
be  exactly  the  same  as  the  corresponding  term  when  NS-EW— 1  (the 
"clockwise"  terms),  except  for  a  change  in  sign  of  H(2,l),  H(l,2) 
and  H(2,3).  The  same  can  be  said  for  the  "counterclockwise"  terms 
with  NS«1,  EW— 1  as  compared  to  those  with  NS«-1,  EW»1. 


h(1,1>ns,ev.i  = 

H(1,1)NS,EV=-1 

H(2,1)NS,EW=1  = 

-H(2,1)NstEW=-l 

h(1,2)ns,ew-i  “ 

'H(1,2)nS,EW— 1 

(134) 

h(2,2)ns,ew-i  = 

H(2,2)NS,EV=-1 

H(1,3)NS,EW-1  * 

h<1,3)ns,ew— 1 

H<2,3)NS,EWol  * 

-H(2,3)ns,ew*_1 

H(1,1)NS=-1,EW=1 

*  H(1,1)NS=1,EW=-1 

h(2,1)ns=-i,ew*i 

-  -H(2,1)MS-1,ew-_1 

h(1,2)ns=-i,ew=i 

=  -H(1,2)ns=1jEW=_1 

(135) 

h(2,2)ns=-i,ew-i 

=  h(2»2>ns=1,EW=-1 

H(1,3)NS— 1,EW-1 

*  H(1,3)NS-1,EW— 1 

H(2,3)NS«-1,EW-1 

"  — H( 2 , 3  >NS=i , EW=— i 

This  reduces  by 

half  (to  12)  the  number  of 

different 

functional  forms  needed  to  describe  this  model.  Furthermore,  from 
the  expressions  (129)  it  can  be  seen  that  generally 


H(l,2)  -  -H(2,3) 


(136) 


Thus,  for  a  given  latitude  and  satellite  direction  of  travel,  the 
velocity  error  sensitivity  of  a  Transit  fix  can  be  described  by 
three  functions:  H(l,2),  H(l,3)  and  H(2,2).  This  reduces  the 
number  of  independent  functions  (of  a  and  X)  to  10  (four  for 
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height  sensitivity  and  six  for  velocity  sensitivity). 

At  first  glance  of  the  numerical  results  shown  in  Section 

9.2  below,  it  appears  that  changing  the  direction  of  travel  (from 
clockwise  to  counterclockwise)  also  only  changes  the  sign  of 
H(2,l),  H(2,2)  and  H(l,3).  Closer  inspection  reveals  that  this  is 
not  exactly  true  however,  and  the  difference  arises  from  equation 
(C15),  which  is  not  homogeneous  in  NS'EV.  The  difference  however 
is  not  significant,  even  at  extreme  latitudes,  as  seen  in  Section 

9.3  below.  Therefore  we  have 


"^clockwise 

o/i  3) 

v  9  'counterclockwise 

l)ciockwise 

-Hf2  1) 

v  9  'counterclockwise 

(137) 

H<2,2)clockwise 

o(2  2) 

'  9  'counterclockwise 

This  reduces 
functions  (again  this 

3  for  height  error 

the  number  of  substantially  different 
is  only  if  the  sign  is  ignored)  to  7,  with 
sensitivity  and  4  for  velocity  error 

sensitivity.  Representative  functions  are: 


clockwise 
H(2,1)  clockwise 

H(*»2)  ciockwise 

H(2,2)  clockwise 
clockwise 
counterclockwise 
H^’2^  counterclockwise 

9.2  Elevation  Angle  Dependence 


(138) 


Figures  14  to  16  illustrate  the  six  components  of  the  H 
matrix  (the  Transit  velocity  and  height  error  sensitivity  model) 
as  functions  of  the  maximum  satellite  elevation  angle  <j,  for  a 
latitude  of  45°  N,  and  a  satellite  direction  of  travel  NS=EV*1 
(satellite  rising  in  the  north  and  passing  east  of  the  receiver). 
Figures  17  to  19  illustrate  intermediate  geometric  quantities  of 
interest,  namely  the  latitude  of  the  satellite  subpoint  Xp,  the 
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bearing  +  from  the  satellite  subpoint  to  the  receiver  and  the 
bearing  ij*  from  the  receiver  to  the  subpoint  (all  at  the  point  of 
closest  approach). 

Figures  20  to  37  give  the  same  set  of  plots  for  the  other 
three  satellite  directions  of  travel.  From  these  plots  it  can  be 
seen  that  the  primary  effect  of  changing  the  direction  of  travel 
is  simply  to  change  the  sign  of  the  error,  as  discussed  in 
Section  9.1  above. 

Furthermore,  from  Figures  14  to  37,  it  can  be  seen  that  at  this 
latitude  the  effect  of  changing  the  satellite  direction  of  travel 
is  to  change  the  sign  of,  but  not  appreciably  change  the  form  of, 
the  dominant  terms  H(l,3)  and  H(2,2).  Only  the  smaller  term 
H(l,2)  significantly  changes  its  form  with  different  directions 
of  travel,  having  a  "clockwise"  form  (with  NE,EV  =  1,1  or  -1,-1) 
and  a  "counterclockwise"  form  (with  NS,EV  -  -1,1  or  1,-1). 
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Figure  18.  Bearing  From  Subpoint  to  Receiver  (NS-EW-1 , lat-45) 


Elevation  Angle 


Figure  19.  Bearing  From  Receiver  to  Subpoint  (NS»EW»l , lat-45) 


Figure  30.  Bearing  From  Subpoint  to  Receiver  (NS--1.EW-1, 1  at-45) 


Figure  31.  Bearing  From  Receiver  to  Subpoint  (NS=-1 , EW=1 , 1 at=45) 
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Figure  36.  Bearing  From  Subpoint  to  Receiver  (NS-EW»-1, lat-45) 


9.3  Latitude  Dependence 


Figures  38  to  69  illustrate  the  dependence  of  the  H 
matrix  elements  (along  with  f  and  ij/  )  on  the  receiver  latitude. 
This  is  done  by  showing  families  of  curves  where  each  curve 
represents  the  elevation  angle  dependence  of  the  variable  of 
interest,  at  a  fixed  latitude.  The  discrete  latitudes  chosen  for 
the  different  curves  of  each  set  are  0°,  20°,  40°,  60°,  70°,  80° 
and  85®  (north  and  south  for  Figures  38  to  45,  north  only  for 
Figures  46  to  69).  Notice  that  the  latitude  spacing  used  is 
greater  at  low  latitudes.  This  is  because  the  effect  of  changing 
latitude  is  much  less  at  low  latitude,  and  in  fact  for  some 
figures  the  0°  and  20°  curves  cannot  be  clearly  separated. 

The  most  notable  difference  at  high  latitudes  is  that  the 
lower  satellite  elevation  passes  do  not  exist.  The  reason  for 
this  is  that  the  polar  point  in  the  Transit  orbit  places  a  lower 
limit  on  the  maximum  elevation  angle  a  (all  Transit  satellites 
are  in  polar  orbits).  This  limit  on  a  is  greater  than  zero  for 
all  receiver  latitudes  above  59°  (as  shown  in  Appendix  A). 

From  Figures  38  to  69  it  can  be  confirmed  that  the  H(2,2) 
function  has  essentially  the  same  form  (except  for  a  change  of 
sign)  for  all  directions  of  travel  at  all  latitudes,  and  that  the 
same  can  be  said  of  H(l,3).  It  can  also  be  verified  that,  as 
predicted  in  Section  9.1  above,  the  terms  H(1 , 2)=-H(2 , 3)  have 
only  two  different  forms  (except  for  a  change  of  sign)  for  all 
directions  of  travel  at  all  latitudes. 

When  negative  latitudes  are  examined  an  additional 
symmetry  appears.  Although  it  is  only  illustrated  here 
explicitly  for  Figures  38  to  45  (the  NS=EW=1  case)  the  general 
symmetry  is  easy  to  imagine  from  the  remaining  Figures  46  to  69, 
and  can  be  proven  using  the  equations  of  Section  8.0.  That  is: 
the  clockwise  families  of  curves  at  negative  latitudes  are  of 
exactly  the  same  form  (ie.  is  the  same  except  for  possibly  a 
change  of  sign)  as  the  corresponding  counterclockwise  families 
for  positive  latitudes: 


(X)clockwise 


±H(i,j) 


(-X) 


clockwise 


(139) 


for  all  i,j.  We  can  see  this  by  following  through  the  equations 
of  Section  8.0,  where  simultaneously  changing  the  sign  of  X  and 
NS*EW  simply: 
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1- changes  the  sign  of  X  in  (C15)  (because  of  NS*EW  on 

the  denominator)  p 

2- leaves  V  unchanged  in  (5) 

3- changes  the  sign  of  cosi^  in  (D6)  (because  X  and  X  both 

change  sign)  p 

4- either  changes  the  sign  of  4  in  (B13)  or  leaves  it 
unchanged  (depending  on  whether  or  not  EV  changes  sign) 

5- either  changes  the  sign  of  sin^  in  (D3)  or  leaves  it 
unchanged  (depending  on  whether  or  not  $  changes  sign) 


In  the  final  equation  (129)  for  the  components  of  H,  the 
net  effect  is  to  change: 


COSl|f 

-*  -COSiJ/ 

sinif; 

-»  -SEWsimJ; 

(140) 

NS'EW 

-►  -NS-EW 

where  SEW  is  negative  if  EW  changes  sign.  The  effect  is  therefore 
to  change  the  sign  of  H(l,l),  H(l,3)  and  H(2,2),  and  to  change 
the  sign  of  the  three  remaining  terms  only  if  EW  changes  sign.  In 
other  words 


H(lf 1) 

^X^ clockwise 

-  ‘"^clockwise 

H(l,3) 

^  ^clockwise 

"  ^clockwise 

(141) 

H(2,2) 

clockwise 

=  -H(2,2)  <-x>clockwise 

H(2, 1 ) 

(X)NS«EW«1  = 

3 

-H(2,l)  (-Mns=1eu=_1 

H(2,l)  <“x)us— l.EW-1 

(142) 

H(2, 1) 

(X)NS-EW— 1  = 

H(2,l) 

— H( 2 » 1)  (-X>NS— 1,EW=1 
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H<1,2)  (X) 


NS-EV-1 

H<1’2>  <X>NS.EV~1 

H(2,3)  (X)ns.EVib1 
H(2,3)  (^ns.ev,.! 


-H(l,2) 

(_X)NS=1,EV=-1 

H(l,2) 

(_X)NS— l,Etf=l 

H(l,2) 

(_X)NS=1,EW=-1 

-H(l,2) 

<-X)NS=-l,EW=l 

i 

to 

fo 

u> 

(-X)NS=l,EW— 1 

H(2,3) 

(-X)NS=-l,EV=l 

H(2,3) 

(-X)NS=l,EV=-l 

-H(2,3) 

(_X)NS— 1,EW-1 

(143) 


(144) 


This  exact  symmetry  nov  replaces  the  approximate  symmetry 
of  equation  (137),  and  leaves  only  five  "basic"  families  (for 
latitudes  from  -90°  to  +90°)  of  functions  of  o,  rather  than  the  7 
families  (some  of  vhich  were  only  approximately  equal,  for 
latitudes  0  to  +90)  listed  in  (138).  These  five  families  can  be 
represented  by: 


clockwise 
0(2 ,X)  clockwise 

0<1’2>  clockwise  <145> 

H(2,2'  clockwise 
H(1,3)  clockwise 


which  are  fully  illustrated  in  Figures  38  to  43. 
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Figure  41.  Latitude  Dependence 
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Figure  50.  Latitude  Dependence  of  H(l,3)  (NS»1,  EW— 1) 


Figure  51.  Latitude  Dependence  of  H(2,3)  (NS*l,  EW=-1) 
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Figure  54.  Latitude  Dependence  of  (NS--1,  EW-l) 


Figure  55.  Latitude  Dependence  of  H(2,l)  (NS=-1,  EW*l) 
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Figure  58.  Latitude  Dependence  of  H(l,3)  (NS=-1,  EW-l) 


Figure  59.  Latitude  Dependence  of  H(2,3)  (NS=-1,  EW=l) 
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Figure  60.  Latitude  Dependence  of  PHI  (NS«-1,  EW«l) 


Figure  61.  Latitude  Dependence  of  PSI  (NS=-l,  Ew=i) 
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10.0  EXPERIMENTAL  VERIFICATION 


The  sensitivity  of  the  Transit  position  fix  to  the  error 
of  the  input  height,  given  by  H(l,l)  and  H(2,l),  is  fairly 
straightforward  and  well  understood,  and  hence  does  not  require 
detailed  verification. 

The  velocity  sensitivities  on  the  other  hand  are  much 
more  complicated  and  hence  warrant  careful  experimental 
verification.  Although  reference  [1]  shows  qualitatively  how  the 
magnitude  of  these  errors  vary  with  elevation  angle  o  at  low 
latitudes,  nothing  could  be  found  in  the  open  literature  to 
indicate  the  sign  of  the  errors,  how  they  vary  with  satellite 
direction  of  travel  (NS,EV),  or  how  they  vary  with  receiver 
latitude.  Therefore  experiments  were  conducted  at  DREO  (latitude 
45°)  in  1983  and  more  extensively  in  1985,  to  measure  these 
errors.  The  1985  results  are  presented  here. 

In  this  1985  experiment  two  stationary  receivers  were 
used  (an  MX1107  and  a  JMR-4)  with  a  -1.0  m/s  velocity 
"measurement"  fed  into  them,  alternately  in  the  north  and  east 
directions.  This  creates  a  velocity  error  (true  -  estimate)  of 
+1.0  m/s  in  the  north  or  east  direction,  to  generate  the 
appropriate  velocity-error-induced  position  errors.  These 
resulting  position  errors  were  recorded,  and  plotted  as  functions 
of  a,  NS  and  EV.  Although  both  receivers  are  dual  channel  (and 
hence  able  to  remove  most  of  the  ionospheric  induced  error), 
there  were  of  course  other  sources  of  position  error  beyond  our 
control  (as  described  in  references  [1]  and  [3]).  These  other 
errors  appear  mostly  as  additive  noise  on  the  error  curves, 
although  these  other  sources  of  error  (especially  the  atmospheric 
induced  sources)  are  not  entirely  uncorrelated  with  the  elevation 
angle.  To  some  extent  these  can  mask  the  details  of  the  smaller 
velocity  induced  terms  (the  clockwise  H(l,2)  and  H(2,3)  terms). 
Therefore  we  use  this  data  here  primarily  to  verify  the  general 
form  of  the  curves  (especially  the  sign).  The  dominant  terms: 
H(l,3)  and  H(2,2)  are  naturally  of  most  concern. 

Figures  70  to  85  present  the  1985  experimental  results. 
Each  figure  shows  three  curves  representing  one  component  of  the 
velocity  error  sensitivity  H  matrix:  the  modelled  version  (solid 
line),  the  JMR  measured  error  (dashed  line)  and  the  MX1107 
measured  error  (dot-dashed  line).  As  expected,  all  errors  clearly 
agree  in  sign  and  general  form  with  the  model.  Even  the  smaller 
terms  can  be  seen  to  agree  on  average  with  the  model's 
prediction,  however  in  this  case  the  noise  is  almost  at  the  same 
level  as  the  deterministic  error. 

It  should  be  mentioned  that  an  adjustment  had  to  be  made 
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to  the  raw  JMR  H(l,2)  and  H(2,3)  data,  in  an  attempt  to  remove  an 
unwanted  error  in  the  recorded  position  data  which  was  not  part 
of  the  actual  position  fix  error.  This  is  an  error  which  is 
introduced  by  the  JMR-4  receiver  after  the  position  fix  is 
performed,  but  before  the  fix  data  is  recorded  for  output.  The 
MX1107  did  not  introduce  this  error  during  our  experiment  because 
it  uses  a  slightly  different  procedure,  as  explained  below. 

The  position  fix  output  from  the  MX1107  is  "at  the  fix 
mark”,  meaning  that  it  is  the  measured  position  of  the  receiver 
at  the  point  in  time  designated  by  the  receiver  as  the  "fix 
mark".  The  fix  mark  is  generally  near  the  centre  of  the  satellite 
pass,  near  the  point  of  closest  approach.  Of  course  this  data  is 
not  available  from  the  receiver  until  after  the  satellite  pass 
has  finished,  and  the  calculations  have  been  performed,  which  is 
generally  about  5  to  9  minutes  after  the  fix  mark.  The  JMR-4,  on 
the  other  hand,  provides  a  position  estimate  that  is  deadreckoned 
from  the  fix  mark  to  the  time  at  which  the  data  is  output.  Thus 
the  JMR-4  position  measurements  have  an  additional  error  in  the 
north  or  east  direction,  of  about  300  to  540  metres,  due  to 
deadreckoning  for  6  to  9  minutes  with  the  1.0  m/sec  velocity 
error.  This  therefore  changed  the  north  position  error  due  to  a 
north  velocity  error,  H(l,2),  and  the  east  position  error  due  to 
an  east  velocity  error,  H(2,3),  which  are  the  small  terms.  In  an 
attempt  to  remove  this  effect,  400  metres  was  added  to  these  two 
terms  (before  plotting  figures  70  to  85).  When  this  error  is 
taken  into  account,  these  terms  are  in  good  agreement  with  both 
the  MX1107  data  and  the  model. 
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Figure  70.  Experimental  Results:  H(l,2)  (NS=EW=1 , lat=45) 


Figure  71.  Experimental  Results:  H(2,2)  (NS=EW=1 , 1at=45) 
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Figure  72.  Experimental  Results:  H(l,3)  (NS=EW»1 , lat«45) 
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Figure  73.  Experimental  Results:  H ( 2 , 3 )  (NS=EW=1 , 1at=45) 


Figure  74.  Experimental  Results:  H ( 1 , 2 )  (NS=1 ,EW=-1 , lat=45) 
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Figure  75.  Experimental  Results:  H(2,2)  (NS=1 ,EW=-i , lat=45) 


Figure  76.  Experimental  Results:  H(l,3)  (NS*1 , EW--1 , lat«45) 


Figure  82.  Experimental  Results:  H(l,2)  (NS=EW=-l ,  lat=45) 


Figure  83.  Experimental  Results:  H(2,2)  (NS=EW=-1 , lat=45) 


11.0  APPLICATIONS 


This  error  model  can  be  used  to  predict  the  accuracy  of 
Transit  fixes  taken  on  a  moving  platform,  as  long  as  the  accuracy 
of  the  velocity  and  height  provided  to  the  receiver  is  known. 

Another  application,  and  the  one  for  which  this  error 
model  was  developed,  is  the  design  of  Kalman  filter  based 
optimally  integrated  multi-sensor  navigation  systems  which  use 
Transit  measurements.  The  H  matrix  described  in  this  report  is 
actually  a  submatrix  of  the  Kalman  filter  "measurement  matrix" 
for  both  the  Marine  Integrated  Navigation  System  (MINS), 
described  in  reference  12 J,  and  the  Primary  Land  Arctic 
Navigation  System  (PLANS),  described  in  [3],  both  of  which  were 
developed  at  DREO. 

Since  the  Transit  orbital  parameters  (H,Vs)  appear 
explicitly  in  the  model,  they  can  be  easily  varied  to  study  the 
effect  of  different  orbits.  This  offers  two  benefits:  firstly  it 
helps  in  selecting  the  optimal  orbit  from  a  system  design  point 
of  view,  and  secondly  it  allows  other  similar  satellite  Doppler 
positioning  systems  to  be  modelled. 

In  the  first  case,  when  choosing  the  best  orbit,  one  must 
account  for  velocity  induced  errors  since  these  do  account  for  a 
large  portion  of  the  total  error  budget.  From  the  equations  of 
Section  8.0,  and  in  particular  equation  (129),  we  can  see  that 
the  most  direct  effect  of  changing  orbital  parameters  comes  from 
the  p/V  term.  This  will  increase  significantly  as  orbital  height 
H  increases,  since  this  directly  increases  the  slant  range  p, 
while  at  the  same  time  decreasing  the  orbital  velocity  Vs  (which 
of  course  is  not  independent  of  H)  which  decreases  V.  Therefore, 
as  one  would  expect,  lower  faster  orbits  lead  to  less  velocity 
error  sensitivity.  However,  a  lower  orbit  also  leads  to  greater 
atmospheric  drag,  which  causes  other  system  problems,  so  to  do 
the  trade-off  analysis  one  must  have  a  quantitative  model,  such 
as  provided  here. 

In  the  second  case,  this  model  can  be  applied  directly  to 
any  similar  Doppler  positioning  system,  such  as  SARSAT  (Search 
and  Rescue  Satellite)  for  example.  SARSAT  works  on  the  same 
Doppler  positioning  principle  as  Transit  does,  and  it  also  has 
circular  orbits,  which  are  almost  polar  (inclination  of  98°). 
This  extra  8°  inclination  can  be  accounted  for  by  the  appropriate 
coordinate  transformation. 

This  error  model  could  also  have  some  more  general  use  as 
a  generic  Doppler  positioning  error  model,  in  that  the  same 
derivations  could  be  followed,  in  the  coordinate  system 
appropriate  to  the  problem. 


-Ill- 


12.0  CONCLUSIONS 


A  detailed  and  general  velocity  error  and  height  error 
sensitivity  model  has  been  derived  in  full  for  dynamic  Transit 
position  fixing,  explicitly  showing  the  position  error  which 
would  result  from  any  error  in  the  velocity  and  height  (which 
must  be  provided  as  inputs  to  the  Transit  receiver  during  the 
position  fix). 

This  model  is  presented  in  algorithmic  form,  and 
numerical  results  are  graphically  presented  to  illustrate  the 
dependence  on  satellite  maximum  elevation  angle,  satellite 
direction  of  travel,  and  receiver  latitude.  Substantial 
experimental  results  are  also  presented,  which  verify  the 
accuracy  of  the  model  for  middle  latitudes  (for  all  satellite 
maximum  elevation  angles  and  all  satellite  directions  of  travel). 

One  somewhat  surprising  result  of  this  general  model  is 
that  at  high  latitudes  the  Transit  error  sensitivity  matrix  is 
significantly  different  than  at  low  or  medium  latitudes.  The 
terms  which  are  small  at  low  latitudes  become  much  larger  (and 
sometimes  change  sign!)  at  high  latitudes,  while  the  normally 
dominant  terms  become  smaller.  Thus  the  error  model  as  previously 
described  (qualitatively)  in  the  literature  was  totally 
inadequate  for  use  at  high  latitudes,  a  deficiency  which  had  to 
be  remedied  for  use  in  the  DREO  developed  Primary  Land  Arctic 
Navigation  System  (PLANS). 
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APPENDIX  A.  Angular  Separation  and  Slant  Range:  6(9),  p( a) 


This  appendix  shows  how  the  angular  separation  0,  between 
the  receiver  and  the  satellite  subpoint,  can  be  expressed  as  a 
function  of  the  satellite  elevation  angle  a.  In  the  process  we 
will  also  obtain  an  expression  for  the  slant  range  p.  This  slant 
range  is  used  in  the  equation  (61)  for  the  measurement  matrix, 
whereas  9  is  needed  to  evaluate  V  and  in  (61). 

Figure  A1  illustrates  the  quantities  of  interest  in  the 
plane  containing  the  receiver,  the  satellite,  the  subpoint  and 
the  earth's  centre.  Applying  the  law  of  cosines  we  obtain: 


p2  «  R2  +  (R+H)2  -  2R(R+H)cos0  (Al) 

and  applying  it  again  we  have: 

(R+H)2  =  p2  +  R1  -  2pRcos(  o+n/2) 

=  p2  +  R2  +  2pRsin<r  (A2) 

sine  =  (H+R)2  -  p2  -  R2  ]  (A3) 

Using  (Al)  in  (A3)  we  have: 

sine  =  {  (H+R)2  -  [  R2  +  (R+H)2  -  2R(R+H)cos0  ]  -  R2  }/2pR 
=  [  2R(R+H)cos0  -  2R2  ]/2pR 

-  [  (R+H)cos0  -  R  ]/p  (A4) 

This  gives  an  expression  for  the  slant  range  p  as  a  function  of 
the  elevation  angle  a  and  separation  angle  9: 
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Figure  Al.  Separation  Angle  in  Vertical  Plane  of  Receiver  &  Satellite. 
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LEGEND: 

0  -separation  angle 
p  -  slant  range 
a  -  elevation  angle 
R  -  earth  radius 
H  -  satellite  height 
r  -  receiver  position 
s  -  satellite  position 
p  -  subpoint  position 


Earth  Center 


5- 


p 


[  (R+H)cos0  -  R  ]/sinc 


(A5) 


This  can  nov  be  squared  and  substituted  in  (Al)  to  obtain  a 
quadratic  equation  for  cos@.  Squaring  (A5)  gives: 

p2  =  [  (R+H)2cos20  +  R2  -  2R(R+H)cos0  ]/sin2o  (A6) 

Substituting  this  in  (Al)  produces: 

(R+H) 2cos2 0  +  R2  -  2R(R+H)cos0 

-  sin2 a[  R2  +  (R+H)2  -  2R(R+H)cos0  1  (A7) 

or 

(R+H)2cos20  -  2R(R+H)cos2 ocos0  +  [R2cos2or  -  (R+H)2sin2o]  *  0 

(A8) 


This  is  quadratic  in  (R+H)cos0,  and  therefore  has  solution: 


(R+H) cos 0  -  2Rcos2c  ±  /(2Rcos2a)2  -  4[R2cos2 a  -(R+H) 2sin2 <x] 

2 


■  Rcos 2  a  ±  /  R2cos2  cr(cos2  a  -1)  +  (R+H)2sin2<r 


«  Rcos2ff  ±  sina  /  (R+H)2  -  R2cos2a  (A9) 


Since  the  satellite  must  be  between  the  horizon  and  the  zenith, 
ve  know  that  0  <  a  <  it/2  and  hence  sina  >  0.  This  positive 
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elevation  angle  also  implies  (as  can  see  from  figure  Al)  that 


(R+H)cos9  >  R 


(A10) 


and  so  the  positive  sign  must  be  chosen  in  (A9).  Therefore: 


cos0  *  R  cos2 q  +  sina  /  1  -  7  R  Vcos^q 
R+H  Ir+hJ 


(All) 


and 


-1 

0  -  cos  (  R  cos2  a  +  sin  a  /  1  - 
l  R+H 

'  R  1 

r+hJ 

2  cos 2  a  J 

(A12) 


The  validity  of  equation  (A12)  can  be  easily  checked  for 
various  limiting  cases.  For  example  if  the  satellite  orbital 
height  H  goes  to  zero  then  the  separation  angle  from  the  receiver 
to  the  satellite  subpoint  0  (at  the  point  of  closest  approach) 
should  approach  zero,  vhich  it  does  according  to  (A12).  Also  if 
the  elevation  angle  a  goes  to  n/2  then  the  separation  angle  0 
should  go  to  zero,  to  vhich  (A12)  does  comply. 

Equation  (A12)  incidentally  places  a  strict  upper  bound 
on  the  separation  angle  0,  at  zero  elevation  angle  a: 


0(a)  <  0(0)  -  cos'1  “  31 


(A13) 


More  importantly,  ve  can  use  (A12)  to  check  the 
consistency  of  the  inputs  a  and  X.  Since  all  Transit  satellites 
are  in  polar  orbits,  all  orbits  contain  the  points  directly  above 
the  poles  (at  orbital  height).  Therefore  the  point  of  closest 
approach  (PCA)  cannot  be  further  than  these  north  or  south  polar 
points.  (If  the  receiver  is  located  exactly  at  a  pole,  then  the 
PCA  vill  be  the  polar  point.)  Therefore  0  cannot  be  greater  than 
the  angle  of  separation  from  the  nearest  pole,  vhich  is  n/2  - 
| X | .  Thus  ve  have 


e  <  |  -  |x| 


(A14) 


so  that 


cos0  >  cos  ^  j  -  |X|)  =  sin|X| 


or 


cos9  >  sin  | X | 


(A15) 


This  is  a  physical  constraint  that,  through  (A12)  effectively 
places  a  lover  limit  on  the  elevation  angle  a  at  high  latitudes. 
Another  interpretation  of  this  limit  is  that  when  the  receiver  is 
close  enough  to  the  pole  so  that  the  satellite  orbital  position 
over  the  pole  is  visible  above  the  horizon,  then  the  maximum 
elevation  angle  a  must  be  at  least  as  great  as  the  elevation  to 
this  polar  orbital  point.  In  view  of  (A13)  we  can  see  that  this 
polar  orbital  point  can  be  seen  anywhere  within  31°  of  the  pole, 
so  that  (A15)  limits  the  value  of  a  for  all  latitudes  above  59° 
(and  below  -59°). 

Now  (All)  can  be  used  in  (A5)  to  obtain  an  expression  for 
the  slant  range  p  directly  as  a  function  of  the  elevation  angle 

<TJ 


p  -  (  Rcos2a  +  sino/(R+H)2  -  R2cos2a  -  RJ/sina  (A16) 


■  R(cos2<r  -  l)/sin<r  +  /R2 (1-cos 2 c)  +  2RH  +  H2  (A17) 


p  -  ✓k2sin2  a  +  2RH  +  H2  -  Rsinc 


(A18) 


This  is  easily  verified  for  a  *  n/2,  in  which  case  the 
satellite  is  directly  overhead  at  range  p-H,  in  agreement  with 
(A18).  It  can  also  be  verified  for  o*0,  in  which  case  the 
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satellite  is  on  the  horizon.  In  this  case  the  range  p  forms  one 
edge  of  a  right  angle  triangle  with  hypotenuse  (R+H)  and  opposite 
side  R.  This  makes 


p  =  /  (R+H) 2  -  R2 

-  /  2RH  +  H2 


(A19) 


which  is  also  in  agreement  with  (A18).  This  coincidentally  places 
bounds  on  the  slant  range  to  the  satellite  at  clossest  approach: 


1,075,000  <  p  <  3,853,716 


(A20) 
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APPENDIX  B.  Heading  of  &  Bearing  from  the  Subpoint: 

«NS,y,  KNS,EV,y 

This  appendix  shows  how  the  latitude  of  the  satellite 
subpoint  X  ,  and  the  heading  of  the  subpoint  path  |3  can  be 
expressed  as  functions  of  the  knovn  quantities  X,  NS  and  EV  (the 
receiver  latitude,  the  satellite  north-to-south  direction  of 
travel  and  the  satellite  east-or-west  direction  from  the 
receiver)  and  the  separation  angle  0.  This  angle  0  can  be 
expressed  in  terms  of  the  knovn  satellite  elevation  angle  a  (as 
shown  in  Appendix  A) . 

The  heading  of  the  subpoint  path  fi  can  be  easily  found  by 
finding  the  north  and  east  components  of  the  subpoint  velocity  in 
the  locally  horizontal  earth  fixed  frame.  We  know  that  the 
satellite  is  in  a  polar  orbit  of  radius  (R+H)  with  orbital  speed 
V  .  Therefore  the  subpoint  moves  directly  "north"  (or  south  if 
N§-1)  in  an  earth  centred  inertial  (non-rotating)  frame  in  an 
"orbit"  of  radius  R  with  speed  RVg/(R+H).  In  the  locally 
horizontal  earth  fixed  frame  the  northerly  component  of  the 
subpoint  velocity  is  therefore  -NS  RVs/(R+H)  and  easterly 
component  of  the  subpoint  velocity  is  due  entirely  to  earth  rate, 
which  is  -RcocosX.  .  Therefore  as  seen  from  Figure  Bl,  the  subpoint 
heading  satisfies: 


tang 


R  cocos  Xp 

NS-R-V  /(R+H) 
s 


(Bl) 


Careful  consideration  of  satellite  direction  of  travel  can  be 
used  to  resolve  the  ambiguity  of  the  inverse  tangent  function, 
leading  to  the  following  solution.  This  consideration  indicates 
that  the  subpath  is  slightly  to  the  west  of  north  if  the 
satellite  is  moving  north  (NS--1)  and  slightly  to  the  west  of 
south  if  the  satellite  is  moving  south  (NS=1),  and  hence: 


(B2) 


This  can  now  be  used  to  express  the  bearing  $  from  the 
subpoint  to  the  receiver  at  the  point  of  closest  approach,  since 
this  bearing  line  (great  circle  from  subpoint  to  receiver)  is 
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orthogonal  to  the  subpoint  path  (a  fact  easily  proven  by 
definition  of  closest  approach)  as  shovn  in  Figure  Bl.  Thus: 


♦  -  0  +  y 


(B3) 


where  direction  of  travel  consideration  (in  Figure  Bl  ve  have 
NS— 1  and  EV-1)  gives: 


V  -  NS-EV  n/2 


(B4) 


therefore: 


♦  =  NS  tan-1  p^n^cosXp)  +  |(1  +  NS  +  NS-EV) 


(B5) 


It  will  turn  out  that  we  actually  need  sin+,  which  is  simpler: 
sin*  -  sinjNS  tan-1  ^R^H^cosXpj  +  £(1  +  NS  +  NS-EV) j  (B6) 


=  cos 


jNS  tan-1  p^KosXp)  +  NS  +  EVnj  j  (B7) 


cos  ^  tan-1  p~2“cos\pj  +  |  +  EVn  j 
-sin  ^  tan-1  ^^^cosX^j  +  EVn  j 


-EVcos  (  tan 


(B8) 


(B9) 


(BIO) 
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This  can  be  further  simplified  by  using  a  trigonometric  identity. 
By  setting 


D  ■  (R+H)  cocos  X  (Bll) 

Vs  P 

ve  have 

cos(tan“*D)  -  1  (B12) 

71+d1 


so  that 


sin$ 


-EV 

/I  +  D2 


(B13) 


J 
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APPENDIX  C.  Latitude  of  the  Subpoint:  Ap(0, A,NS,EW) 


This 

subpoint  X„ 


appendix  shows  how 
can  be  expressed 


the  latitude  of  the  satellite 
as  a  function  of  the  known 
quantities  ‘'X,  NS  and  EV.  To  do  this  we  consider  the  spherical 
triangle  shown  in  Figure  Cl,  with  vertices  at  the  receiver,  the 
satellite  subpoint  (at  the  point  of  closest  approach)  and  the 
north  pole.  Ue  can  consider  X  and  6  to  be  known,  and  we  seek  to 
find  X  .  (In  Appendix  D  we  shall  also  find  an  expression  for  \f/. ) 
From  tHe  law  of  cosines  for  spherical  triangles,  where  the 
"sides"  as  shown  in  Figure  Cl  are  treated  as  angles  subtended  at 
the  earth's  centre,  we  have: 


cos(n/2-A)  =  cos9cos( it/2-A  )  +  sin6sin( n/2-X  )cos(2n-t) 

P  P  (Cl) 


sinX  =  cosSsinXp  +  sin0cosXpcost  (C2) 


Now  from  equation  (B5)  we  can  expand  cost,  and  using  simple 
trigonometric  identities,  simplify  it: 


cost  =  cos  (  NS  tan  1  r(R+H)«o  _v  +  it  +  NSii(l+EV)  ^  (C3) 

l  lnr-cosXPJ  2  2  J 


cost  =  -sin  ^  NS  tan-1  ^(R+tOoo.^  j  +  NS(n  +  EVn)  j  (C4) 


>  -NSsin  f  tan  1r(R+H)w  .  )  +  it  +  EVit  'k  (C5) 

(  l-V^cosXpJ  2  2  J 


-  -NScos  ^  tan-1  |p4jO«cosX  j  +  EVjt  j  (C6) 
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NS-EVsinf  tan"1  ((R+H)«  x  ^  ^ 

(  nr— cosXpJ  J 


Now  the  argument  of  tan"  in  (C7)  is  a  small  angle  (the  angle 
between  the  satellite  subpoint  path  and  the  meridian).  Its 
maximum  value  is  at  the  equator,  (where  the  horizontal  earth  rate 
is  maximum),  where  it  is  approximately  4  degrees.  Therefore  to  a 
good  approximation  we  can  write: 


cos*  «  NS-EVf  (R+H)«  „  'l 

l  v“c0spj 


This  can  now  be  used  in  (C2)  to  create  a  quadratic  equation  in 

sinX  : 

P 


sinX  =  cosGsinX  +  NS*EW(R+H)wsin9cos2  X 


cos  0s  in Xp  +  NS*EW-Ksin0(l-sin2  Xp)  (CIO) 


where 


K  ■  (R+H)«/Vs  =  0.074416 


(Cll) 


Thus  (CIO)  can  be  written: 


(NS-EW-Ksin0)sinIXp  -  cosGsinX^  +  (sinX  -  NS-EV-KsinG)  =  0 


which  is  quadratic  in  sinXp,  and  hence  has  solutions: 


sinX,  »  cos8  ±  /cos29  -  4NS*EW>Ksin9(sinX  -  NS’EW’KsinG) 
p  2NS-EW-Ksin0 
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(C13) 


Choosing  the  positive  or  negative  sign  is  not  difficult.  As  was 
seen  from  equation  (A13),  the  separation  angle  0  cannot  be 
greater  than  32  degrees,  so  that  cos@  >  0.84  and  sin6  <  0.53.  Ve 
also  sav  in  (Cll)  above  that  K  is  small  (<0.075).  This  implies 
that 

cos9  >  11  (C14) 
2Ksin6 


Therefore  the  negative  sign  in  equation  (C13)  yields  the  only 
possible  solution.  Also  as  the  separation  angle  0  becomes  small, 
ve  must  have  X  approaching  X,  as  is  the  case  with  the  negative 
sign.  Thus  the  latellite  subpoint  latitude  X  is  (note  that  NS2  = 
EW2  *  1  and  that  by  definition  X  c  [-jt/2,  n/ZJ )  s 


X  =  sin  *  fcosO  -  /cos 2 8  -  4Ksin8(NS»EVsinA  -  Ksin0) 
P  l  2NS-EW-Ksin0 


where  K  is  given  by  (Cll). 


(C15) 


Nov  ve  can  find  a  simple  approximation  to  equation  (C15) 
by  ignoring  higher  order  terms  in  K: 


sinX  a  cos8  -  /cos28  -  4K-NS,EVsin8sinX)  (C16) 

P  2NS-EW-Ksin0 


Using  the  Taylor  series: 


(x+e)%  -  x%  +  e-fcx"%  +  e2 * [->Ax~3/2]  +  ...  (C17) 

2 


ve  have: 


/cos 2 0  -  4K-NS-EWsin0sinX  =  cos0  -  4K-NS-EWsin8sinX  (C18) 

2cos0 


-2K*NS*EWtan9sinX 


(C19) 
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Thus  (C16)  becomes: 


sinX  *  2K*NS*EWtanQsinX 

p  2NS-EW-Ksine 


(C20) 


sinX  (C21) 

COS0 


This  approximation  is  included  here  for  its  intuitive  value,  and 
to  provide  a  simple  test  of  reasonableness  which  can  be  applied 
to  the  solution  of  (C15). 
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APPENDIX  D.  Bearing  to  Subpoint: 


*<e,Xp,\) 


i  Now  the  spherical  triangle  in  Figure  Cl  can  be  used  to 

find  the  other  angles.  This  Figure  shows  the  situation  in  which 
the  satellite  is  east  of  the  receiver  (EV=1).  For  the  other  case 
.  (EV*-1)  the  Figure  is  simply  reflected  about  the  receiver 

meridian.  In  this  case  -♦  becomes  4  and  the  bearing  angle  t|>  from 
north  is  negative.  Therefore  the  law  of  sines  for  spherical 
triangles  can  be  applied  to  solve  for  the  bearing  to  the 
receiver  at  the  subpoint.: 


so 


sin(rt/2  -  X) 
EW-sin(-4) 


sin(tt/2  -  Xp) 
EW*  sintj/ 


cosX  *  cosXp 

-sin+  sin<j/ 


(Dl) 


(D2) 


» 


sin  ^  =  -sinfcosXp 
cosX 


(D3) 


Applying  the  law  of  cosines  for  sides  gives: 


cos(n/2-Xp)  =  cos9cos(  Jt/2-X)  +  sin9sin(  re/2-X)cos^  (D4) 


or 


sinX 

P 


cosQsinX  +  sinQcosXcos^ 


(D5) 


so  that 


J 


cos^  *  sinXp  -  cos8sinX 
sinQcosX 


(D6) 
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